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Abstract. Spectral analysis of hot luminous stars requires ad- 
equate model atmospheres which take into account the effects 
of NLTE and radiation driven winds properly. Here we present 
significant improvements of our approach in constructing de- 
tailed atmospheric models and synthetic spectra for hot lumi- 
nous stars. Moreover, as we regard our solution method in its 
present stage already as a standard procedure, we make our pro- 
gram package WM-basic available to the community (down- 
load is possible from the URL given below). 

The most important model improvements towards a realis- 
tic description of stationary wind models concern: 
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(i) A sophisticated and consistent description of line block- 
ing and blanketing. Our solution concept to this problem 
renders the line blocking influence on the ionizing fluxes 
emerging from the atmospheres of hot stars — mainly 
the spectral ranges of the EUV and the UV are affected 
— in identical quality as the synthetic high resolution 
spectra representing the observable region. In addition, 
the line blanketing effect is properly accounted for in the 
energy balance. 

(ii) The atomic data archive which has been improved and 
enhanced considerably, providing the basis for a detailed 
multilevel NLTE treatment of the metal ions (from C 
to Zn) and an adequate representation of line blocking 
and the radiative line acceleration. 

(iii) A revised inclusion of EUV and X-ray radiation pro- 
duced by cooling zones which originate from the sim- 
ulation of shock heated matter. 

This new tool not only provides an easy to use method for 
O-star diagnostics, whereby physical constraints on the prop- 
erties of stellar winds, stellar parameters, and abundances can 
be obtained via a comparison of observed and synthetic spec- 
tra, but also allows the astrophysically important information 
about the ionizing fluxes of hot stars to be determined automat- 
ically. Results illustrating this are discussed by means of a basic 
model grid calculated for O-stars with solar metalUcity. To fur- 
ther demonstrate the astrophysical potential of our new method 
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we provide a first detailed spectral diagnostic determination of 
the stellar parameters, the wind parameters, and the abundances 
by an exemplary application to one of our grid-stars, the 09.5Ia 
0-supergiant a Cam. 

Key words: Line: formation - Stars: atmospheres - Stars: 
early type - Stars: mass-loss - Stars: individual: a Cam - X- 
rays: stars 



1. Introduction 

Spectral analyses of hot luminous stars are of growing astro- 
physical interest as they provide a unique tool for the determin- 
ination of the properties of young populations in galaxies. This 
objective, however, requires spectral observation of individual 
objects in distant galaxies. That this is feasible has already been 
shown by Steidel et al. (1996) who detected galaxies at high 
redshifts (z ^ 3.5) and found that the corresponding optical 
spectra show the typical features usually found in the UV spec- 
tra of hot stars. In order to determine stellar abundances and 
physical properties of the most UV-luminous stars in at least the 
Local Group galaxies via quantitative UV spectroscopy another 
principal difficulty needs to be overcome: the diagnostic tools 
and techniques must be provided. This requires the construc- 
tion of detailed atmospheric models and synthetic spectra for 
hot luminous stars. It is a continuing effort of several groups to 
develop a standard code for solving this problem. (Basic papers 
of the different groups are Hillier and Miller (1998), Schaerer 
and de Koter (1997) and Pauldrach et al. (1994, 1994a, 1998).) 

The most important output of this kind of model calcula- 
tion are the ionizing fluxes and synthetic spectra emitted by the 
atmospheres of hot stars. As these spectra consist of hundreds 
of not only strong, but also weak wind-contaminated spectral 
lines which form the basis of a quantitative analysis, and as 
the energy distribution from hot stars is also used as input for 
the analysis of emission line spectra (e. g., of gaseous nebulae) 
which depend sensitively on the structure of the emergent stel- 
lar flux, a sophisticated and well tested method is required to 
produce these data sets accurately. However, this is not an easy 
task, since modelling hot star atmospheres involves replicating 
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a tightly interwoven mesh of physical processes: the equations 
of radiation hydrodynamics including the energy equation, the 
rate equations for all important ions (from H to Zn) including 
the atomic physics, and the radiative transfer equation at all 
transition frequencies have to be solved simultaneously. 

The most complicating effect in this system is the over- 
lap of thousands of spectral lines of different ions. Especially 
concerning this latter point we have made significant progress 
in developing a fast numerical method which accounts for the 
blocking and blanketing influence of all metal lines in the en- 
tire sub- and supersonically expanding atmosphere. As we have 
found from previous model calculations that the behaviour 
of most of the spectral lines depends critically on a detailed 
and consistent description of line blocking and line blanketing, 
(cf. Pauldrach et al., 1994; Sellmaier et al., 1996; Taresch et 
al., 1997; Haser et al., 1997), special emphasis has been given 
to the correct treatment of the Doppler-shifted line radiation 
transport, the corresponding coupling with the radiative rates 
in the rate equations, and the energy conservation. 

In Section ^ we will demonstrate that the realistic and 
consistent description of line blocking and blanketing and the 
involved modifications to the models lead to changes in the 
energy distributions, ionizing continua, and line spectra with 
much better agreement with the observed spectra when com- 
pared to previous, not completely consistent models. This will 
obviously have important repercussions for quantitative analy- 
sis of hot star spectra. 

In the next two sections we will first summarize the general 
concept of our procedure and then discuss the current status of 
our treatment of hydrodynamical expanding atmospheres. 

2. The general method 

The basis of our approach in constructing detailed atmo- 
spheric models for hot luminous stars is the concept of ho- 
mogeneous, stationary, and spherically symmetric radiation 
driven winds, where the expansion of the atmosphere is due 
to scattering and absorption of Doppler-shifted metal lines 
(Lucy and Solomon (1970~). In contrast to previous papers of 
this series (cf. Pauldrach et al., (1994, 1994a, 1998) and papers 
referenced therein) the above approximations are now the most 
significant ones for the present approach — others of similar 
importance have meanwhile been dropped (see below). These 
approximations are, however, quite restrictive, since only the 
time-averaged mean of the observed spectral features can be 
described correctly by our method. Nevertheless we are confi- 
dent that it is reasonable to continue with the stationary, spher- 
ically symmetric approach and to improve its inherent physics 
since the detailed comparison with the observations, which is 
the only way to demonstrate the reliability of this concept, leads 
to promising results (cf. Section 4). 

Before we describe the latest improvements in detail we 
first summarize the principal features of our procedure of sim- 
ulating the atmospheres of hot stars. (For particular points a 
comprehensive discussion is also found in the papers cited 
above.) 



Fig. |l] gives an overview of the physics to be treated in vari- 
ous iteration cycles. A complete model atmosphere calculation 
consists of three main blocks, 

(i) the solution of the hydrodynamics 

(ii) the solution of the NLTE-model (calculation of the radi- 
ation field and the occupation numbers) 

(iii) the computation of the synthetic spectrum 

which interact with each other. 

In the first step the hydrodynamics is solved in dependence 
of the stellar parameters (effective temperature Toff, surface 
gravity log 5, stellar radius i?* (defined at a Rosseland opti- 
cal depth of 2/3), and abundances Z (in units of the corre- 
sponding solar values)) and of pre-specified force multiplier 
parameters (k, a., 5), which are used for describing the radia- 
tive line acceleration. In addition, the continuum force is ap- 
proximated by the Thomson force, and a constant temperature 
structure (T(r) = Tcff) is assumed in this step. In a second step 
the hydrodynamics is solved by iterating the complete contin- 
uum force gcii") (which includes the opacities of all important 
ions) and the temperature structure (both are calculated using 
a spherical grey model), and the density p{r) and the velocity 
structure v{r) in a pre-iteration cycle. In a final outer iteration 
cycle these structures are iterated again together with the line 
force obtained from the spherical NLTE model. (New 

force multiplier parameters, which are depth dependent if re- 
quired, are deduced from this calculation.) [] 

The main part of the code consists of the solution of the 
NLTE-model. In this step the radiation field (represented by 
the Eddington-flux Hi,{r) and the mean intesity J^{r)), the fi- 
nal temperature structure, occupation numbers rii (r), and opac- 
ities Xiy and emissivities rji, are computed using detailed atomic 
models for all important ions. For the solution of the radia- 
tive transfer equation the influence of the spectral lines (i. e., 
the UV and EUV line blocking) is properly taken into account 
in addition to the usual consideration of continuum opacities 
and source functions which consist of Thomson-scattering and 
free-free and bound-free contributions of all important ions. 
Moreover, the shock source functions produced by radiative 
cooling zones which originate from a revised simulation of 
shock heated matter are included also. For the calculation of 
the final NLTE tem.perature structure the line blanketing effect, 
which is a direct consequence of line blocking, is considered 
by demanding luminosity conservation and the balance of mi- 
croscopic heating and cooling rates. The rate equations which 
yield the occupation numbers contain collisional (C^ ) and ra- 
diative (i?ij) transition rates, as well as low-temperature dielec- 
tronic recombination and Auger-ionization due to K-shell ab- 
sorption (considered for C, N, O, Ne, Mg, Si, and S) of soft 
X-ray radiation arising from shock-heated matter. (Further de- 
tails concerning the solution method of the NLTE-model are 
described in Section ||.) 

The last step consists of the computation of the synthetic 
spectrum. In dependence of the occupation numbers, the opac- 

' This latter step is currently not available for the download version 
of the code; it will be made available for version 2.0. 
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Fig. 1. Schematic sketch of a complete model atmosphere calculation. Starting procedures are presented in brackets. For a discussion see the 
text. 
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ities and the emissivities, a complete synthetic spectrum is 
computed using a formal integral solution of the transfer equa- 
tion in the observer's frame. To accurately account for the vari- 
ation of the line opacities and emissivities due to the Doppler 
shift, the calculation is performed on a properly adapted spatial 
microgrid which effectively resolves individual line profiles. 
All lines which, through their Doppler shift, can contribute 
to a given frequency point are considered, (cf. Puis and Paul- 
drach (1990), Pauldrach et al. (1996)). 

As results of the iterative solution of this system of equa- 
tions we obtain not only the synthetic spectra and ionizing 
fluxes which can be used in order to determine stellar parame- 
ters and abundances, but also the hydrodynamical structure of 
the wind (thus, constraints for the mass-loss rate M and the 
velocity structure v{r) can be derived). 

3. The consistent NLTE model 

The construction of realistic models for expanding atmospheres 
requires a correct and extremely consistent description of the 
main part of the simulation, the NLTE model, which, in addi- 
tion, should not suffer from drastic approximations. From our 
continuing effort to come to a reasonable approach of a solu- 
tion to this problem it turned out that the most crucial point in 
our present treatment is an exact description of line blocking 
and blanketing. 

The effect of line blocking — mainly acting in between the 
He II and the H I edge — is that it influences the ionization and 
excitation and the momentum transfer of the radiation field sig- 
nificantly. This of course has important consequences for both 
the spectral line formation and the dynamics of the expanding 
atmosphere. Nevertheless, it is still not a common procedure 
to treat the line opacities and emissivities in the radiative trans- 
fer equation and their back-reaction on the occupation numbers 
via the radiative rates correctly. We will therefore first discuss 
the effects of line blocking and blanketing for expanding atmo- 
spheres of hot stars in more detail. 

The huge number of metal lines present in hot stars in the 
EUV and UV attenuate the radiation in these frequency ranges 
drastically by radiative absorption and scattering processes (an 
effect known as line blocking). Only a small fraction of the 
radiation is re-emitted and scattered in the outward direction; 
most of the energy is radiated back to the surface of the star pro- 
ducing there a backwarming. Due to the increase of the Rosse- 
land optical depth (tross) resulting from the opacities enhanced 
by the line blocking, and, in consequence, of the temperature, 
the radiation is redistributed to lower energies (this refers to 
line blanketing). In principle these effects influence the NLTE 
model with respect to: 

(i) the radiative photoionization rates Rik, 

(ii) the radiative bound-bound rates Rij , 

(iii) the radiation pressure ^rad, 

(iv) the energy balance. 

The terms of the first two items are directly connected to the 
radiation field, and line blocking in general reduces them con- 



siderably. Concerning the third item, the blocked incident radi- 
ation reduces the radiative acceleration term in the inner part, 
whereas it can be enhanced in the outer part due to multiple 
scattering processes (cf. Puis (1987) and references therein). In 
contrast to this, the energy equation — last item — is mostly 
influenced by the impact of the line opacities, and this blanket- 
ing effect results in an increased temperature (steeper gradient) 
in the deeper layers of the photosphere. 

Although the method for treating blanketing effects is well 
established for cold stars, where the atmospheres are hydro- 
static and where the assumption of LTE is justified (cf. Ku- 
rucz, 1979 and 1992), the work to develop an adequate method 
for hot stars, where not only NLTE effects are prominent, but 
where the atmospheres are also rapidly expanding, is still under 
way. In this — our — case, in addition to the four items given 
above, the solution of the radiative transfer also has to account 
for the lineshift caused by the Doppler effect due to the veloc- 
ity field. The important effect of this point is that the velocity 
field increases the frequential range which can be blocked by a 
single line (see below). In the presence of a velocity field the 
blocking effect is therefore more pronounced. 

Concerning the basic requirements for calculating adequate 
line opacities and source functions for expanding atmospheres 
of hot stars we have to concentrate on the following points: 

(1) consistent NLTE occupation numbers, 

(2) a complete and accurate line list in connection with de- 
tailed atomic models, 

(3) a proper concept for treating the line blocking with due 
regard to the lineshifts in the wind, in the course of which 
the method for solving the complete radiative transfer in- 
cluding the spectral lines has to be efficient with regard to 
computational time, 

(4) a correct treatment of the influence of the blanketing ef- 
fect on the temperature structure, 

(5) an adequate approximation of the EUV and X-ray radia- 
tion produced by cooling zones of shock-heated matter. 

3.1. The concept of the solution of ionization and excitation 

It is obvious that ionization and excitation plays the major role 
in calculating the emergent flux and spectrum of a hot star. 
Therefore, a consistent and accurate description of the occu- 
pation numbers is extremely important for a realistic solution 
of the NLTE model. 

Figure ^ presents a sketch of our iteration scheme for the 
calculation of the occupation numbers. 

In dependence of the abundances (Z), the density {p{r)) 
and velocity {v{r)), and a pre-specified temperature structure 
(Tg{r)) (see section^), the occupation numbers are determined 
by the rate equations containing collisional (Cy ) and radiative 
(Rij) transition rates. The most crucial dependency of the rates 
is not the density, which is nevertheless important for the col- 
lisional rates and the equation of particle conservation, but the 
velocity field which enters not only directly into the radiative 
rates via the Doppler shift, but also indirectly through the ra- 
diation field determined by the equation of transfer, which in 
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Fig. 2. Iteration scheme for the calculation of the NLTE occupation numbers. An Accelerated Lambda Iteration procedure is involved in the 
blocking and blanketing cycles. Note that two successive iteration cycles (I,II) of different quality are applied for the blocking and blanketing 
part of the model calculation (see text). 
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turn is again dependent on the Doppler shifted line opacities 
and emissivities. 

For the calculation of the radiative bound-bound tran- 
sition probabilities Rij we make use of the Sobolev-plus- 
continuum method (Hummer and Rybicki, 1985; Puis and 
Hummer, 1988). Only for some weak second-order lines in the 
subsonic region of the atmospheric layers where the continuum 
is formed might this be just a poor approximation (cf. Sellmaier 
et al., 1993). A more important point of our procedure concerns 
the problem of self-shadowing (cf. Pauldrach et al., 1998). This 
problem occurs because the rate equations are not really solved 
simultaneously with the radiative transfer, but instead in the 
framework of the accelerated lambda iteration (ALT), in which 
the radiation field and the occupation numbers are alternately 
computed (cf. Pauldrach and Herrero, 1988). Hence, the radia- 
tion field which enters into a bound-bound transition probabil- 
ity is already affected by the line itself, since the line has also 
been considered for the blocking opacities. This procedure will 
lead to a systematic error if a line transition dominates within 
a frequency interval (see section 3.3.1 ). The solution for cor- 
rectly calculating the bound-bound rates even in these circum- 
stances is quite simple and has been described by Pauldrach et 
al. (1998, section 3.2). 

The spherical transfer equation yields the radiation field at 
2,500 frequency points (see below) and at every depth point, in- 
cluding the layers where the radiation is thermalized and hence 
the diffusion approximation is a proper boundary condition. 
The solution includes all relevant opacities. In particular, the 
effects of wind and photospheric EUV line blocking on the 
ionization and excitation of levels are treated on the basis of 
4 million lines, with proper consideration of the influence of 
the velocity field on the line opacities and emissivities and on 
the radiative rates. 

Regarding the latter point, the inclusion of line opacities 
and emissivities in the transfer equation, two different concepts 
are employed for iterating the occupation numbers and the tem- 
perature structure until a converged radiation field (r) and 
H^{r) is obtained. In a first step, a pre-iteration cycle with an 
opacity sampling method is used (method I). This procedure 
has the advantage of only moderate computing time require- 
ments, allowing us to perform the major part of the necessary 
iterations with this method. Its disadvantage, however, is that it 
involves a few substantial approximations (cf. section 3.3). In 
a second step, the final iteration cycle is therefore solved with 
the detailed radiative line transfer (method II). Although this 
procedure is extremely time-consuming, it has the advantage 
that it is not affected by any significant approximations. With 
this second method, blocking factors Bj {r, v) and Bh {r, v) 
are calculated, defined as the ratio of the radiative quantities 
obtained by considering the total opacities and emissivities to 
those which include only the corresponding continuum values 
(cf. Pauldrach et al., 1996). Bj{r,v) and Bh{i\i^) are then 
used as multiplying factors to the continuum quantities calcu- 
lated in the next NLTE-ALI-cycle with the current continuum 
opacities, in order to iterate the radiative rates Rij (both con- 



tinuum and lines) and the resulting occupation numbers until 
convergence (details are described in section 3.3). 

In total, almost 1000 ALI iterations are required by the 
complete NLTE procedure, divided into blocks of 30 itera- 
tions each. (One iteration comprises calculation of the occu- 
pation numbers and the radiation field.) Up to 31 of these iter- 
ation blocks are performed using the opacity sampling method 
(method I), updating the temperature structure and the Rosse- 
land optical depth after each third ALI-iteration, and the total 
opacities and emissivities after each iteration block. The fol- 
lowing iterations are then all performed using method II, up- 
dating temperature, optical depth, and opacities and emissiv- 
ities as before, and additionally calculating the blocking fac- 
tors after each iteration block. (Several iteration blocks using 
method II can be executed, but 1 is usually sufficient — see 
below.) In this phase the radiative transfer solved in the ALI- 
iterations is just based on continuum opacities and emissivities, 
and the blocking factors are applied to get the correct radiative 
quantities used for calculating the radiative rates. 

As a final result of the complete iteration cycle, the con- 
verged occupation numbers, the emergent flux, and the final 
NLTE temperature structure are obtained. 

3.2. The atomic models 

It is obvious that the quality of the calculated occupation num- 
bers and of the synthetic spectrum is directly dependent on the 
quality of the input data. We have therefore extensively revised 
and improved the basis of our model calculations, the atomic 
models. 

Up to now the atomic models of all of the important ions 
of the 149 ionization stages of the 26 elements considered 
(H to Zn, apart from Li, Be, B, and Sc) have been replaced 
in order to improve the quality. This has been done using the 
Superstructure program (Eissner et al., 1974; Nussbaumer and 
Storey, 1978), which employs the configuration-interaction ap- 
proximation to determine wave functions and radiative data. 
The improvements include more energy levels (comprising a 
total of about 5,000 observed levels, where the fine structure 
levels have been "packed" together^ and transitions (compris- 
ing more than 30,000 bound-bound transitions for the NLTE 
calculations and more than 4,000,000 lines for the line-force 
and blocking calculations]^'^ and 20,000 individual transition 



" Note that artificial emission lines may occur in the blocking cal- 
culations if the lower levels of a fine structure multiple! are left un- 
packed but the upper levels of the considered lines are packed — the 
fine structure levels of an ionization stage should either be all packed 
or all unpacked. 

Note further that the consistency of the model calculation requires 
that the wavelength of the bound-bound transition connecting packed 
levels in the NLTE calculations to be identical to the wavelength of 
the strongest component of the multiplet considered in the blocking 
calculations in order to solve the line radiative transfer and especially 
the problem of self-shadowing properly. 

■* Our line list also includes transitions to highly excited levels 
above our limit of considering the level structure; occupation numbers 
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Table 1. Summary of Atomic Data. Columns 2 and 3 give the number of levels in packed and unpacked form; in columns 4 and 5 the number 
of lines used in the rate equations and for the line-force and blocking calculations are given, respectively. 

levels lines levels lines 



Ion 


packed 


unpacked 


rate eq. 


blocking 


Ion 


packed 


unpacked 


rate eq. 


blocking 


Cii 


36 


73 


284 


11005 


S VI 


18 


32 


59 


142 


Cm 


50 


90 


520 


4406 


S VII 


14 


26 


39 


1031 


C IV 


27 


48 


103 


229 


Ar V 


40 


86 


328 


3007 


C V 


5 


7 


6 


57 


Ar VI 


42 


93 


400 


1335 


N III 


40 


82 


356 


16458 


Ar VII 


47 


87 


483 


2198 


Niv 


50 


90 


520 


4401 


Ar VIII 


15 


27 


41 


111 


N V 


27 


48 


104 


229 


Mn III 


50 


141 


364 


175593 


N VI 


5 


7 


6 


57 


Mn IV 


50 


124 


467 


131821 


II 


50 


117 


595 


39207 


Mn V 


50 


124 


508 


61790 


Oiii 


50 


102 


554 


24506 


Mn VI 


13 


25 


35 


87 


IV 


44 


90 


435 


17933 


Fe II 


50 


148 


405 


227548 


V 


50 


88 


524 


4336 


Fe III 


50 


126 


246 


199484 


VI 


27 


48 


102 


231 


Fe IV 


45 


126 


253 


172902 


Ne IV 


50 


113 


577 


4470 


Fe V 


50 


124 


451 


124157 


Ne V 


50 


110 


534 


2664 


Fe VI 


50 


138 


452 


60458 


Ne VI 


50 


112 


343 


1912 


Fe VII 


22 


62 


91 


10123 


Mg III 


50 


96 


529 


2457 


Fe viii 




yo 




^ 1 1 i 


Mg IV 


50 


117 


589 


3669 


Co III 


50 


141 


469 


200637 


iVlg V 




1 nn 


JH- / 




V^U 1 V 




Q7 


/U 




Tvlg VI 


21 


44 






Co V 








1 89780 
1 oz. / ou 


Al IV 


50 


96 


529 


2523 


Co VI 


43 


113 


317 


124053 


Al V 


50 


117 


588 


18317 


Co VII 


34 


80 


246 


50270 


Al VI 


19 


37 


41 


153 


Ni III 


40 


102 


281 


131508 


Si III 


50 


90 


480 


4044 


Ni IV 


50 


146 


528 


183267 


Si IV 


25 


45 


90 


245 


Niv 


41 


97 


70 


179921 


Si V 


50 


98 


531 


3096 


Ni VI 


45 


126 


253 


186055 


Si VI 


50 


116 


596 


3889 


Ni VII 


43 


113 


317 


123386 












Ni VIII 


34 


80 


246 


43778 


P V 


25 


45 


90 


245 












P VI 


14 


26 


41 


1096 


Cu IV 


50 


124 


477 


17466 












Cu V 


50 


146 


527 


30457 


S V 


44 


78 


404 


903 


Cu VI 


50 


126 


246 


10849 



probabilities of low-temperature dielectronic recombination 
and autoionization). 

Additional line data were taken from the Kurucz (1992) line 
list and from the Opacity Project (cf. Seaton et al., 1994; Cunto 
and Mendoza, 1992); the latter was also a source for photoion- 
ization cross-sections (almost 2,000 data sets have been incor- 
porated). CoUisional data have become available through the 
IRON project (see Hummer et al., 1993) — almost 1,300 data 
sets have been included. 

Table |l| gives an overview of the ions affected by the im- 
provements. (Users of the program package WM-basic should 
note that the model calculations will become inconsistent if the 
atomic data sets are changed haphazardly by those who are not 
familiar with the source code.) 



of these upper levels are estimated using the two-level approximation 
on the basis of the (known) occupation number of the lower level. 



3.3. The treatment of line blocking 

As the thermal width of a UV metal line covers just a few mA, 
a simple straightforward method would require considering ap- 
proximately 10^ frequency points in order to resolve the lines 
in the spectral range affected by line blocking. Such a proce- 
dure would lead to a severe problem concerning the computa- 
tional time. The alternatives are either to calculate the complete 
radiative transfer in the comoving frame — again a time con- 
suming procedure — , or to use a tricky method which saves a 
lot of computation time through the application of some minor 
approximations (method I), dropping these approximations in 
the final iteration steps (method II) in order to come to a re- 
alistic solution. Our treatment described here uses the second 
approach. 

Although frequently applied, a method using opacity dis- 
tribution functions (ODFs — cf. Labs, 1951; Kurucz, 1979), 
where the opacities are rearranged within a rough set of fre- 
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quency intervals in such a way that a smoothly varying func- 
tion is obtained which conserves the statistical distribution of 
the opacities, is not applicable in our case, since there is no ap- 
propriate way to treat the lineshift in the wind, and due to the 
rearrangement of the opacities the frequential position of the 
lines is changed. This, however, prevents a correct computa- 
tion of the bound-bound transitions used for the solution of the 
statistical equilibrium equations. 

The approach best suited for our purpose in the first step 
(method I) is the opacity sampling technique (cf. Peytre- 
mann, 1974; Sneden et al., 1976) which compared to the ODF- 
method is computationally a bit more costly, but does not suffer 
from the limitations mentioned above. This method allows us 
to account for the lineshift in the wind and the correct influence 
of line blocking on the bound-bound transitions (cf. Section ||, 
item (ii)), since it preserves the exact frequential position of the 
lines. 

3.3.1. The opacity sampling method 

Following the idea of the opacity sampling, a representative set 
of frequency points is distributed in a logarithmic wavelength 
scale over the relevant spectral range, and the radiative trans- 
fer equation is solved for each point. (For O-stars the actual 
range depends on Tcff; for hot objects the lower value is at 
« 90 A and for cooler objects the upper value is at 2000 A; 
note that accurate ionization calculations require extending the 
line blocking calculations to the range shortward of the He II 
edge — cf. Pauldrach et al., 1994) 

In this way the exact solution is reached by increasing 
the number of frequency points. A smooth transition is ob- 
tained when the number of frequency points is increased up 
to the number — 10^ — which is required to resolve the ther- 
mal width of the UV lines. It is obvious however that conver- 
gence can be achieved already with significantly less points 
(see below). Furthermore, special blocking effects on selected 
bound-bound transitions can be investigated more thoroughly 
by spreading additional frequency points around the line tran- 
sition of interest. 

In the following subsection we will investigate how many 
sampling points are required in order to represent the physical 
situation in a correct way. 

The influence of line blocking on the photoionization integrals. 
The most important effect of line blocking on the emergent 
spectrum is the influence on the ionization structure via the 
photoionization integrals 



(1) 



This can be verified from Fig. ^ where it is shown that the mean 
radiation field changes rapidly over the frequency interval 
covered by a typical smooth bound-free cross section tjik — 
several 100 A are affected. (Note that dielectronic resonances 
which may occur in addition are not shown here.) 
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Fig. 3. Mean radiation field Ju together with the photoionization 
cross section ait of the ground state of H (in arbitrary units). 



1.10 



1 .05 




0.90 



100 1000 10000 

number of sampling points 

Fig. 4. Accuracy of the normalized photoionization integral FUk of 
the groundstate of H in dependence of an increasing number of sam- 
pling points within the Lyman continuum. 



It is obvious from Fig. || that the photoionization rates are 
sensitive functions of the blocking influence on J^, and hence, 
on the number of sampling points in the relevant frequency 
range. 

In order to determine the number of sampling points re- 
quired for an accurate description we performed empirical tests 
by evaluating typical photoionization integrals in dependence 
of an increasing numbers of sampling points. Fig. ^ shows as an 
example the normalized photoionization integral of the ground 
state of hydrogen. The result is that 1,000 sampling points 
within the Lyman continuum guarantee a sufficient accuracy 
of about 1 to 2 percent. By means of a separate investigation 
Sellmaier (1996) showed the given number of sampling points 
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to be reasonable, since it reproduces the actual line-strength 
distribution function quite well. 

The treatment of the lineshift. The total opacity at a certain 
sampling frequency v is given by adding the line opacity xib to 
the continuum opacity Xc 



Xv = Xib(i') ^Xc{v), 



(2) 



where xvo is the sum over all (integrated) single line opacities 
Xi multiplied by the line profile function tf\(y) 



lines 



Here X\ is 

X\ = ^ (niBtj - rijBji) , 
and the analogous expression for the emissivity is 



(Br 



1 

An 



(3) 



(4) 



(5) 



Bji, und Aji are the Einstein coefficients of the line tran- 



sition at the frequency i/p, and h is Planck's constant.) 

In the static part of the atmosphere a line's opacity covers 
with its (thermal and microturbulent) Doppler profile (/^d only 
a very small interval around the transition frequency j/q (illus- 
trated in Fig. ^ on the right hand side of both figures; note that 
with regard to our sampling grid about 40 percent of the avail- 
able lines are treated in this part). The effect of these lines on 
the radiation field is nevertheless considerable (cf. Fig. ||), if 
the lines are strong enough to become saturated. 

In the expanding atmospheres of hot stars the effect of line 
blocking is enhanced considerably in the supersonic region due 
to the nonlinear character of the radiative transfer A velocity 
field v{r) enables the line to block the radiation also at other 
frequencies v = vqII + v{r)/c), i.e., the Doppler shift in- 
creases the frequency interval which can be blocked by a single 
line to a factor of w 100. On the other hand, the velocity field 
reduces the spatial area where a photon can be absorbed by a 
line. If a line is optically thick, however, the effect of blocking 
will ultimately be increased compared to a static photosphere. 

The lineshift due to the velocity field is applied to the in- 
dividual line opacities before the summation in eq. || is carried 
out at each sampling and depth point (otherwise the effect of 
the lineshift would be underestimated with respect to the ratio 
of line width to sampling distance — see below). However, in 
our approach this is done by applying the Doppler shift of the 
radial ray to all p-rays (see Fig. ignoring the angular depen- 
dence of the Doppler shift (see below). Apart from the intrinsic 
character of the sampling method this is the most restrictive ap- 
proximation in our first iteration cycle. The frequency shift is 
thus treated in a statistical sense, its main effect — increase of 
the frequential range of line absorption — is nevertheless in- 
cluded properly and that is what has to be iterated in this first 
cycle. 



o 





frequency points 

Fig. 5. upper panel: simply shifting the line profile along vcmf 
(represented by the curve) at each radius grid point (standard opac- 
ity sampling) causes the line to be missed at most frequency points; 
lower panel: this problem is solved by assuming a boxcar profile 
for each depth point with a width corresponding to the difference in 
Doppler shift between to successive radius points ("Doppler-spread 
opacity sampling"). 




to observer 



Fig. 6. (p, 2;)-geometry for the spherically symmetric radiative trans- 
fer. For any given depth point, a different Doppler shift must in prin- 
ciple be considered for every p-ray, since the projected velocity varies 
with cos 9. However, as no analogy to the boxcar method exists in this 
case, in our opacity sampling method we take the Doppler shift of the 
central ray as being representative for all other rays (see text). 



From the upper panel of Fig. ^ it is obvious that if the 
line opacity is simply shifted along the comoving frame fre- 
quency (vcmf) to every radius point successively, many fre- 
quency points will miss the line, since the radius grid is too 
coarse to treat large lineshifts in the observer's frame. This be- 
haviour is corrected by convolving the intrinsic Doppler profile 
of the line with a boxcar profile ip^^y representing the velocity 
range around each radius point (Fig. |], lower panel). 
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The boxcar profile is the mean profile obtained by consider- 
ing the velocity shifts Av of the two corresponding intermesh 
points (vi, V2) on both sides of the regarded radius grid point 
in the way that the gaps in the frequency grid are closed. This 
can be expressed in terms of the Heaviside function 9: 



{v2 — v) — 6{vi — v) 



(6) 



vi and V2 are the observer's frame frequencies belonging to 
the velocities of two successive radius points (ri and r2), i. e., 
1^1,2 = fo(l + w(ri_2)/c). Assuming thermal Doppler broaden- 
ing for the intrinsic line profile. 



with X 



^0 



(7) 



where A^d is the thermal Doppler width, the convolution 
{ifD (g) (^Ai;)(i^) results in the final profile function 



(j){v) = (<^D <8) ^PAv){v) 



crf(a;2 — x) — erf(xi — x) 

2{X2 - Xi)Avy) 



(8) 



This profile can be used for the entire sub- and supersonic re- 
gion. For Aw < wthcrm it gives, as a lower limit, the ordinary 
opacity sampling, and for sufficiently high velocity gradients 
(Aw > wthorm) the integration over a radius interval represents 
the Sobolev optical depth (rsob(?')) of a local resonance zone 
for a radial ray 



Ar = 



Xi 



Xi 



Xi 



l-(-l) 



2{X2 - xi) ■ 
r2 — ri c 



AfD 



{r2 - ri) 



V2 



^0 
Tsoh{r). 



Av 
dr 



(9) 



Note that at sufficiently high velocity gradients all lines are 
included in the radiative transfer if the sampling grid is fine 
enough (see also Sellmaier (1996)). In this case our Doppler- 
spread opacity sampling method therefore becomes an exact 
solution. 

In principle, one would have to account for the angular vari- 
ation of the Doppler shift in a similar manner (see Fig. |^), but 
as no analogy to the boxcar profile method exists in this case, 
as mentioned above we simply apply the (correctly calculated, 
radially dependent) opacities {x\h{v^ f)) of the central p-ray to 
the other p-rays, regarding these opacities as being represen- 
tative. A welcome result of this simplification is the fastness 
of the method, a very important consideration in this iteration 
cycle. 

Concerning our WM-basic program package running on a 
normal scalar processor, however, the method is still not fast 





Fig. 7. Logarithm of the ionization fractions of N III {upper panel), 
N IV {middle panel), and N V (lower panel) versus density and iter- 
ation block number for an O supergiant model (Toff = 29 000 K, 
log.g = 3.0, R^/Rq = 27.0). The region between two successive 
contour lines (one iteration block) corresponds to 30 iterations. 
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enough (a model calculation would require an amount of com- 
puting time of about 20 hours). The reason is that the Rybicki- 
method which is used in this step for the solution of the second- 
order form of the equation of transfer (cf. Mihalas, 1978) re- 
quires more than 80% of the computing time of a model cal- 
culation. (Note that the Rybicki-method is applied in each it- 
eration just once per frequency point; in order to improve the 
accuracy, the radiative quantities are then further iterated in- 
ternally by using the moments equation of transfer (cf. Mi- 
halas, 1978). Because of strong changes in the opacities and 
emissivities within the NLTE iteration cycle it is necessary to 
start with the Rybicki-method nevertheless.) We have therefore 
rethought the solution concept of the Rybicki-scheme and de- 
veloped a method which is 10 times faster on a vector processor 
and 3 to 5 times faster on a scalar processor — the actual factor 
depends on the quality of the level-2 BLAS functions available 
with professional compiler programs and which do most of the 
work in our method (see 

In order to illustrate the behaviour of convergence of our 
method 1, the ionization fractions of N III, IV, and V are shown 
versus density and the iteration block number in Fig. ^ for 
the first 600 iterations as an example. As displayed, the model 
convergences within 400 iterations — the remaining iterations 
are requ ired to warrant the luminosity conservation (see sec- 



tion 



3.5). The steep increase of Nv in the wind part results 



from the EUV and X-r ay radiation produced by shock-heated 
matter (see section 3.6). 



We finally note that first results obtained with a prelimi- 
nary version of this procedure have already been published. 
Sellmaier et al. (1996) showed that their NLTE line-blocked 
O-star wind models solve the longstanding Ne III problem of 
H ll-regions for the first time, and Hummel et al. (1997) carried 
out NLTE line-blocked models for classical novae. 



Special problems. From first test calculations performed in the 
manner described we recognized and solved two additional 
nontrivial problems: 

The first problem concerns the artificial effect of self- 
shadowing (see above) which occurs because the incident in- 
tensity used for the calculation of a bound-bound transition 
that enters into the rate equations is already affected by the line 
transition itself, since the opacity of the line has been used for 
the computation of the radiative quantities in the previous it- 
eration step. If the lines contained in a frequency interval are 
of almost similar strength, this is no problem, since the used 
intensity /^^^ (r) calculated at the sampling point represents a 
fair mean value for the true incident radiation of the individual 
lines in the interval. If however, a line has a strong opacity with 
a dominating influence in the interval, the intensity taken at the 
sampling point for the same bound-bound transition in the ra- 
diative rates is much smaller than the true incident radiation 
for this line, because the line has already influenced this value 
considerably. In consequence the source function of this line 
is underestimated and the radiative processes — the scattering 




Fig. 8. Part of the synthetic EUV spectrum (200-1600 A) of the S-45 
supergiant model (Toff = 45000 K, logg = 3.6, R,/Rq = 18.0) 
calculated with the Doppler-spread opacity sampling method. The up- 
per panel shows the spectrum resulting from the standard treatment of 
the transfer equation which produces several strong artificial emission 
lines. The lower panel shows the correct spectrum resulting from our 
modified treatment. 

part is mostly affected — are not correctly described in the way 
that the line appears systematically too weak. 

The solution to this problem is rather simple: in calculat- 
ing the bound-bound rates of the dominating lines, we use an 
incident intensity which is independent of the lines in the con- 
sidered interval (cf. Pauldrach et al., 1998). 

The second problem involves the discretization of the trans- 
fer equation in its differential form, for computing the radiative 
quantities (Feautrier method). In the standard approach (see, 
for example, Mihalas, 1978) the equation of transfer is written 
as a second-order differential equation with the optical depth r 
as the independent variable: 



(10) 



where S is the source function and u = + )> with /+ 
and being the intensities in positive and negative t direction 
along the ray considered. 
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This differential equation is then converted to a set of dif- 
ference equations, one for each radius point i on the ray. 



dr2 



du I 



i + 



du I 
dr I 



Uj+l 



resulting in a linear equation system 

a^Mi-i + biUi + CiUi+i — Si 



with coefficients 



Ci = - (|(Ti+l - ri_l)(Ti+l - Ti)) 

bi = 1 - a,; - Ci , 



(11) 



(12) 



(13) 



(14) 



(and appropriate boundary conditions). This linear equation 
system has a tridiagonal structure and can be solved econom- 
ically by standard linear-algebra means.^ Note that the equa- 
tions contain only differences in r, which can easily be calcu- 
lated from the opacities and the underlying z-grid (cf. Fig. 0) 



n+i 



(15) 



with Xi being the opacity at depth point i. 

The equation systems are well-behaved if the opacities and 
source functions vary only slowly with z. Caution must be 
taken if this cannot be guaranteed, for example, whenever a 
velocity field is involved at strong ionization edges or with 
the opacity sampling method at strong lines, since the veloc- 
ity field shifts the lines in frequency, causing large variations 
of the opacity from depth point to depth point for a given fre- 
quency. In particular, a problematic condition occurs if a point 
with a larger-than-average source function Si and low opacity 
Xi borders a point with a high opacity Xi+i (™d low or av- 
erage source function Si+i). In reality, this large source func- 
tion should have little impact, since it occurs in a region of low 
opacity, and thus the emissivity is small. However, the structure 
of the equations is such that the emission is computed to be on 
the order of 



A/ w • Ar 

~ \{Si+i + Si) ■ \{xt+i + Xi){zi+i 



(16) 



where, if the other quantities are comparatively small (in ac- 
cordance with our assumptions), the term SiXi+i dominates,^ 

^ In practice, a Rybicki-type scheme (cf. Mihalas, 1978; and Ap- 
pendix |a[ this paper) is used for solving the equation systems for all 
p-rays simult aneo usly, since the source function contains a scattering 
term (see eq. 4.4 ) which redistributes the intensity at each radius shell 
over all rays intersecting that shell. 

* The physical reason for the failure of the system is that the source 
function only has meaning relative to its corresponding opacity. Mul- 
tiplying the source function from one point with the opacity at another 
point is complete nonsense. 



leading to artificially enhanced emission. In Figure ^ (upper 
panel) we show the exaggerated emission of the strongest spec- 
tral lines in the emergent flux of a stellar model computed using 
this standard discretization, leading to false results. Even a sim- 
ple example can serve to illustrate this effect, as demonstrated 
in Appendix B. 

However, with a subtle modification of the equation system 
coefficients the method can nevertheless be salvaged. The sub- 
tle point involves writing the transfer equation as an equation 
not in T, but in z for derivation of the coefficients, since only 
this formulation treats correctly the z-dependence of x- 

l±(i^)=u-S (17) 
Xdz \x<izj 

(Note that the grid should still be spaced so as to cover t more- 
or-less uniformly.) Again approximating the differential equa- 
tion with a system of differences we obtain 



1 

Xi 



d 

dz 



1 du 
X dz 



1 



1 du 

X dz 



f 1 du 



Zi 



1 - z. 



1 



Ui+1 - u 



1 



1 X 



Zi+l ^ Zi 



XiA-l Zi 



+ Zi) — ^{zi + Zi_i) 



so that 

ai = 

Ci = 
b, = 



{\Xi{zi+\ 
{\Xi{zi^\ 



-i) 
-i) 



Xi,i-\i.Zi 

X'i+\Azi+\ 



z^)y 



(18) 



(19) 



(20) 



1 



Even though these coefficients seem not too different from 
those of the standard method, their impact on the computed 
radiation field is significant, as witnessed by the drastic im- 
provement in the emergent flux shown in the lower panel of 
Figure H. The crucial difference in the coefficients is that the 
first factor in a and c now contains only the local opacity. (We 
naturally make the corresponding changes in the coefficients of 
the moments equation as well.) 

Test calculations have shown that for the second factor in 
the coefficients the geometric mean (an arithmetic mean on a 
logarithmic scale) 



Xi+l^^ = • Xi (21) 

gives good results, as demonstrated in Figure^, where the spec- 
trum of a model computed with the opacity sampling method 
is compared to that of our detailed radiative line transfer, de- 
scribed in the next section. Considering the relative coarseness 
of the opacity sampling method, and the fact that the detailed 
line transfer suffers none of the approximations of the sam- 
pling method, the agreement is indeed remarkable. Note again 
that through our single-p-ray approximation for the sampling 
opacities (see above), our method I (opacity sampling) cannot 
produce P Cygni profiles, since the P Cygni emission is a direct 
result of the different Doppler shifts of a particular spectral line 
along different rays. 
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Fig. 9. Comparison of the synthetic EUV spectrum (200-1600 A) of the S-45 supergiant model (Tcff = 45000 K, log g = 3.6, R, /Rq = 
18.0) calculated with the Doppler-spread opacity sampling method (thick line) and the detailed method (thin line) on the same frequency grid. 
Due to simplifications in our implementation the sampling method cannot produce P Cygni emission; nevertheless it provides an extremely 
good basis for the final iterations using the detailed method. 



3.4. The detailed radiative line transfer 

The detailed radiative line transfer (method II) removes the 
two most significant simplifications of our opacity sampling 
method (method I), i. e., it accounts for: 

( 1 ) Correct treatment of the angular variation of the opacities, 

(2) Spatially resolved line profiles^] (implying coiTect treat- 
ment of multi-line effects). 

Whereas in method I the former is completely ignored, the 
lack of spatial resolution was already compensated for to a 
large extent through the use of our Doppler-spread sampling. 
(Multi-line interaction is partly included in our method I, but 
without regard for the sign of the Doppler shift (using just 
that of the central ray), and without regard for the order of 

' Note that this will not by itself solve the problem of self-shad- 
owing, since that is an intrinsic property of any method using an "in- 
cident radiation" in solving for the bound-bound radiative rates with a 
continuum already affected by the transition being considered. In the 
iteration cycle using method II we therefore also have to apply our 
correction for self-shadowing. 



the lines along the ray within a radius interval, as the Doppler- 
spread sampling effectively "maps" the lines to the nearest ra- 
dius point.) 

With all major approximations removed, the biggest short- 
coming that remains in method II is that only Doppler broad- 
ening is considered for the lines, as Stark broadening has not 
yet been implemented. However, this is of no relevance for the 
UV spectra, as it concerns only a few lines of Hydrogen and 
Helium in the optical frequency range. It will, however, be im- 
portant for our future planned analysis of the optical H and 
He lines. (Stark broadening is not considered in the sampling 
method either, but here this is of minor significance, as all other 
approximations are much more serious.) Method II would re- 
main a sampling in frequency if the frequential resolution were 
chosen so coarse that the whole Doppler profile of a line would 
fit between two frequency points; only if the velocity field car- 
ries a spectral line over a point in the frequency grid will that 
line be considered for radiative transfer (at the coiTect radial 
position). Nevertheless (with the exception of missing Stark 
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Fig. 10. Diagram (not to scale) illustrating the basic relationship of the rest-frame frequencies of spectral lines (i^cMp) to observer's frame 
frequency (yobs) for one particular (non-core) p-ray in the spherically symmetric geometry (cf. Figure |^. Shown are two spectral lines which 
get shifted across the observer's frame frequency by the velocity field in the wind. The dots represent the stepping points of the adaptive 
microgrid used in solving the transfer equation in the detailed radiative line transfer. 



broadening), at all points of the frequency grid the radiative 
transfer is solved for correctly. 

In contrast to our method I, where the symmetry and our 
assumption of only radially (not angular) dependent Doppler 
shifts allowed solving the transfer equation for only one quad- 
rant]^ a correct treatment of the both red and blue Doppler- 
shifted line opacities (see Figure ||) requires a solution in two 
quadrants^ (corresponding to, from the observer's viewpoint, 
the front and back hemispheres; the rotational symmetry along 
the line-of-sight is taken care of through the angular integration 
weights). 

The method employed is an adaptation of the one described 
by Puis and Pauldrach (1990), using an integral formulation of 
the transfer equation and an adaptive stepping technique which 
ensures that the optical depth in each step ("microgrid") does 
not exceed Ar = 0.3, so that the radiation transfer in each 
micro-interval can be approximated to high accuracy by an an- 
alytical formula assuming a linear run of opacity and emissivity 
between the micro-interval endpoints: 

/(ro) = r S{T)e-(^-^°^ At + I{T,,)e~^^--^''\ (22) 

J To 

where the integral is performed as a weighted sum on the mi- 
crogrid 



E 



-ir.-ro) 



(23) 



The 2nd-order differential representation of the transfer equation 
accounts for both the left- and right-propagating radiation simultane- 
ously, the unknowns being the symmetric averages of the two. 

' A one-quadrant solution is also possible, but requires both a red- 
and a blue-shifted opacity for each {jp, 2)-point and separate treatment 
of the left- and right-directed radiation, thus being equivalent in com- 
putational effort to the two-quadrant solution that solves for radiation 
going in only one direction. 



each "microintegral" being evaluated as 

" ^(T)e-(^-^') dr = w'f^Sin) + wf^S{n+^) (24) 
with weights 



(b) 



.-At 



-At, 



(25) 



where Ar^ — r^+i — r^. Note that in constructing the opaci- 
ties and emissivities, all line profile functions ipi (cf. eq. |^ are 
evaluated correctly for the current microgrid-(z,p)-coordinate 
on the ray. Only the slowly-varying occupation numbers (or 
equivalently, the integrated, frequency-independent line opac- 
ities Xi and emissivities !]{) and the velocity field are interpo- 
lated between the regular radius grid points. 

Figure |l^ depicts schematically the relationship between 
the Doppler-shifted frequencies of spectral lines (which are 
constant in the comoving frame) and the observer's frame fre- 
quency for which the radiative transfer is being calculated. The 
figure also illustrates the line overlap in accelerating, expand- 
ing atmospheres: lines clearly separated in the comoving frame 
(slices parallel to the (z^, x)-plane) overlap in the observer's 
frame (slice parallel to the (2, x)-plane at i^obs) due to large 
Doppler shifts many times the intrinsic (thermal and microtur- 
bulent) linewidth. The areas shaded in dark gray correspond to 
the spatially resolved Sobolev resonance zones of the two lines 
for this particular observer's frame frequency and p-ray. Note 
that the dimensions are not to scale, i. e., the intrinsic width of 
the lines, and consequently the thickness in z of the resonance 
zones, has been greatly exaggerated in relation to the total ve- 
locity shift. 

All lines whose maximum Doppler shift Ai' = ±vqVoo/c 
puts them in range of the observer's frame frequency for which 
the radiative transfer is being calculated are considered for that 
frequency point. In Figure |T^, these correspond to those lines 
whose rest frequencies lie in the gray band in the {y, z) -plane 
at 2: = 0. 
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Fig. 11. Flux of our S-29 supergiant model (Toff = 29 000 K, log g = 3.0, R, /Rq = 27.0) after 1 and 5 iteration blocks of method II. 



After the occupation numbers have converged in the itera- 
tion cycle using method I, one iteration block with method II is 
usually sufficient for full convergence of the model, as demon- 
strated by Figure |ll], where the emergent spectrum of our S-30 
model after 1 iteration block of method II is compared to the 
spectrum resulting from 5 iteration blocks. 

A high-resolution spectrum is computed for the purpose of 
comparison with observations (wavelength range usually from 
900 to 1600 A) after full convergence of the model. This spec- 
trum is generated with exactly the same procedure as used for 
the detailed line blocking calculations, and would thus cor- 
respond to a further iteration block of method II. The high- 
resolution spectrum is then merged with the lower-resolution 
blocking flux (A < 900 A) for the final flux output (Figure |l2|). 

3.5. Line blanketing 

Line absorption and emission also has an important effect on 
the atmospheric temperature structure. The corresponding in- 
fluence on the radiation balance is usually referred to as line 
blanketing. The objective now is to calculate an atmospheric 
temperature stratification which conserves the radiative flux 
and which treats the impact of the line opacities and emissiv- 



ities properly. In principle there are three methods for calcu- 
lating electron temperatures in model atmospheres. The com- 
monly used one is based on the condition of radiative equilib- 
rium. The second one uses a flux correction procedure, and the 
third one is based on the thermal balance of heating and cool- 
ing rates. As the first method has some disadvantages (see be- 
low), we use the second and the third method (cf. Pauldrach et 
al., 1998). In deeper layers (tross > 0.1) where true absorptive 
processes dominate we use the flux correction procedure, and 
the thermal balance is used in the outer part of the expanding 
atmosphere (tross < 0.1), where scattering processes start to 
dominate. 



The flux correction procedure. The idea of this method is 
straightforward: The local temperature has to be adjusted in 
such a way that the radiative flux is conserved. This requires, 
however, that the temperature is the dominant parameter on 
which the flux depends, and that the effect of a change in tem- 
perature on the flux is known. The first condition is certainly 
the case for tross > 0.1. With regard to the second condi- 
tion, a law for the temperature structure is required which is 
controlled by some global parameters that can be adjusted in a 
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Fig. 12. The final — completely converged — spectrum from 200 to 2300 A of the S-45 supergiant model. The observable UV region redwards 
of the Lyman edge has been computed with a high resolution (for comparison with observed spectra). 



proper way in order to conserve the flux. The "Hopf function", 
usually applied for the grey case, has an appropriate functional 
dependence which has been adapted to the spherical NLTE case 
by Santolaya-Rey et al. (1997) recently in a general way. The 
main characteristic of the Hopf method is that the Rosseland 
optical depth is the decisive parameter on which the tempera- 
ture stratification depends. Thus, in deeper layers the temper- 
ature structure can be calculated efficiently by using this new 
concept of the NLTE Hop/function: 



4 TRo 



(?n(tRoss) + TRoss) 



(26) 



where tross is the radial optical depth in the spherical case, 



dfRoss = XRoss(?-) 

and 

g'N(rRoss) ^ 900 + {go - Qoo) exp(-7fRoss) 



(27) 



(28) 



is the spherical NLTE Hopf function, where the parameters 

Qq, Qoc, and 7 are fitted to a predefined run of the (/nCtross) 
stratification (cf. Santolaya-Rey et al., 1997). Test calculations 



performed with fixed parameters (qo, q^o, and 7) and without 
metal lines lead to almost identical results for the temperature 
structures obtained by our code and the completely indepen- 
dently developed code of Santolaya-Rey et al. (cf. Pauldrach et 
al., 1998). The reliability of the method has been further proven 
by the resulting flux conservation which turns out to be on the 
1% level. 

In the next step, line blocking has to be treated consistently. 
Although the line processes involved are complex, they always 
increase the Rosseland optical depth (tross)- In the deeper lay- 
ers (tross > 0.1) this leads directly to an enhancement of 
the temperature law (backwarming). Using the method of the 
NLTE Hopf functions we thus have to increase the parameter 
go first by using the flux deviation at tross ~ 3. This param- 
eter is updated in the corresponding iteration cycle until the 
flux is conserved at this depth on the 0.5% level. Afterwards 
the same is done with the parameter q^o at an optical depth 
of tross ~ 0.1. In case the flux deviation at tross ~ 3 be- 
comes larger than 0.5% the parameter go is iterated again with 
a higher priority. As a last step in this procedure the parameter 
7 is adjusted in order to conserve the flux at an optical depth of 

"TRoss ~ 1- 



A.W. A. Pauldrach, T. L. Hoffmann, M. Lennon: NLTE line blocking and blanketing 



17 



,x1 O 



6^1 



2 



4- 



A 



O 



A 
O V 



4 




Fig. 13. The temperature structure versus the Rosseland optical depth and the iteration block number for the same O supergiant model as in 
Fig. |(Tcff = 29 000 K, logp = 3.0, ii./ii© = 27.0). 



The resulting temperature structure and the corresponding 
flux deviation for an O supergiant model (Totf — 29 000 K, 
log g = 3.0, i?*/i?0 = 27.0) ai-e shown in Fig. |l| and Fig. 0, 
respectively. As can be inferred from the figure, the flux is con- 
served for this model with an accuracy of a few percent. (We 
note that from test calculations where line blocking was treated, 
but parameters of the NLTE Hopf function were held fixed, thus 
effectively ignoring blanketing effects, we found a flux devia- 
tion which already starts in the inner part (tross < 50) and 
reaches a value of up to 50% at tross ~ 0.1.) This clearly 
shows the importance of blanketing and backwarming effects 
and the need to include them. As these test calculations have 
also shown that absorptive line opacities dominate the total 
opacity down to an optical depth of tross > 0.1, the tempera- 
ture structure is influenced by backwarming effects in the entire 
atmosphere — cf. eq. 

The thermal balance. In the outer part of the expanding atmo- 
sphere (tross < 0.1), where scattering processes start to domi- 
nate, the effects of the line influence on the temperature struc- 
ture are more difficult to treat. Of the two possible treatments, 
calculating for radiative equilibrium or for thermal balance, we 



have chosen the latter one, as the convergence of the radia- 
tive equilibrium method turned out to be problematic since the 
tross -values are small in this part, and, hence, most frequency 
ranges are optically thin. (This has recently also been proven 
by Kubat et al. (1999), where the corresponding equations of 
the method are also presented.) In calculating the heating and 
cooling rates (Hummer and Seaton, 1963), all processes that 
affect the electron temperature have to be included — bound- 
free transitions (ionization and recombination), free-free tran- 
sitions, and inelastic collisions with ions. For the required iter- 
ative procedure we make use of a linearized Newton-Raphson 
method to extrapolate a temperature that balances the heating 
and cooling rates. 



Fig. 13 displays for the model above the resulting temper- 
ature structure vs. the number of iterations and shows a pro- 
nounced bump and a successive decrease of the temperature 
in the outer atmospheric part. (Note that the mismatch of the 
heating and cooling rates which immediately goes to 0% in the 
outer part (tross < 0.07) where it is applied for correcting the 
temperature structure has already been presented by Pauldrach 
et al. (1998).) 
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3.6. Revised inclusion of EUV and X-ray radiation 

The EUV and X-ray radiation produced by cooling zones 
which originate from the simulation of shock heated matter 
arising from the non-stationary, unstable behaviour of radia- 
tion driven winds (see Lucy and Solomon (1970), who found 
that radiation driven winds are inherently unstable, and Lucy 
and White (1980) and Lucy (1982), who explained the X-rays 
by radiative losses of post-shock regions where the shocks are 
pushed by the non-stationary features) is, together with K-shell 
absorption, included in our radiative transfer. The primary ef- 
fect of the EUV and X-ray radiation is its influence on the ion- 
ization equilibrium with regard to high ionization stages like 
N V and O VI (cf. the problem of "superionization", the detec- 
tion of the resonance lines of O VI, N V, S VI in stellar wind 
spectra (cf. Snow and Morton, 1976); in a first step this problem 
was investigated theoretically by Cassinelli and Olson, 1979) 
where the contribution of enhanced direct photoionization due 
to the EUV shock radiation is as important as the effects of 
Auger-ionization caused by the soft X-ray radiation (cf. Paul- 
drach, 1987; Pauldrach et al., 1994 and 1994a). In order to treat 
this mechanism accurately it is obviously important to describe 
the radiation from the shock instabilities in the stellar wind flow 



properly. Note that in most cases a small fraction of this radia- 
tion leaves the stellar wind to be observed as soft X-rays with 
Lx/L\,o\ ~ 10"'^ (cf. Chlebowski et al., 1989). Thus, the reh- 
ability of the shock description can be further demonstrated by 
a comparison to X-ray observations, by ROSAT for instance. 

In principle, a correct calculation of the creation and de- 
velopment of the shocks is required for the solution of the 
problem. This means that a detailed theoretical investigation 
of time-dependent radiation hydrodynamics has to be per- 
formed (for exemplary calculations see Owocki, Castor, and 
Rybicki (1988) and Feldmeier (1995)). However, these calcu- 
lations favour the picture of a stationary "cool wind" with em- 
bedded randomly distributed shocks where the shock distance 
is much larger than the shock cooling length in the acceler- 
ating part of the wind. They also indicate that only a small 
amount of high velocity material appears with a filling factor 
not much larger than / « 10^^, and jump velocities of about 
Us = 300 . . . 700 km/s which give immediate post-shock tem- 
peratures of approximately Tg = 1 x 10^ to 8 x 10^ K. We 
also note that the reliability of these results was already demon- 
strated by a comparison to ROSAT-observations (cf. Feldmeier 
et al., 1997). 
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On the basis of these results we had developed an empiri- 
cal approximative description of the EUV and X-ray radiation, 
where the shock emission coefficient 

4(0 = -fnpn,A,{Ts{r) n,) (29) 

was incorporated in dependence of the volume emission coef- 
ficient Au calculated by using the Raymond and Smith (1977) 
code for the X-ray plasma, the velocity-dependent post shock 
temperatures T^, and the filling factor / which enter as fit pa- 
rameters — these values are determined from a comparison of 
the calculated and observed ROSAT "spectrum". With this de- 
scription the effects on the high ionization stages (N V, O VI) 
lead to synthetic spectral Unes which reproduce the observa- 
tions almost perfectly (cf. Pauldrach et al., 1994 and 1994a). 
However, with this method we were not able to reproduce the 
ROSAT-observations with the same model parameters simulta- 
neously (see below). We therefore had to determine the filling 
factor and the post-shock temperatures by a separate and hence 
in view of our concept not consistent procedure (cf. Hillier et 
al., 1993). In order to overcome this problem refinements to our 
method are obviously required. 

In the present treatment the outlined approximative descrip- 
tion of the EUV and X-ray radiation has been revised. The ma- 
jor improvement consists of the consideration of cooling zones 
of the randomly distributed shocks embedded in the stationary 
cool component of the wind. Up to now we had assumed, for 
reasons of simplicity, that the shock emission is mostly char- 
acterized by the immediate post-shock temperature, i. e., we 
considered non-stratified, isothermal shocks. This, however, 
neglects the fact that shocks have a cooling structure with a 
certain range of temperatures that contribute to the EUV and 
X-ray spectrum. Our revision comprises two modifications to 
the shock structure. The first one concerns the inner region of 
the wind, where the cooling time can be regarded to be shorter 
than the flow time. Here the approximation of radiative shocks 
can be applied for the cooling process (cf. Chevalier and Ima- 
mura, 1982). The second one concerns the outer region, where 
the stationary terminal velocity is reached, the radiative accel- 
eration is negligible, and the flow time is therefore large. Here 
radiative cooling of the shocks is of minor importance and the 
cooling process can be approximated by adiabatic expansion 
(cf. Simon and Axford (1966), who investgated a pair of re- 
verse and forward shocks that propagate through an ambient 
medium under these circumstances). For our purpose we fol- 
lowed directly the modified concept of isothermal wind shocks 
presented recently by Feldmeier et al. (1997). 

Compared to eq. ^ we account for the density and tem- 
perature stratification in the shock cooling layer by replacing 
the values of the volume emission coefficient (K,j{Ts{r) ric)) 
through adequate integrals over the cooling zones denoted by 
A^{Ts{r)). Thus, e^(r) is replaced by 

etW = 7-«pneA,(T;(r)), (30) 



where 

r^Xs 

A,(T,(r))=±l / P{r')AAnir')-gir'))dr\ (31) 

Xs J 

r 

and r is the location of the shock front, r' is the cooling length 
coordinate with a maximum value of Xg, the plus sign corre- 
sponds to forward and the minus sign to reverse shocks, and 
f{r') and g{r') denote the normalized density and temperature 
structures with respect to the shock front. The improvement of 
our treatment is now obviously directly connected to the de- 
scription of the latter functions. In the present step we used the 
analytical approximations presented by Feldmeier et al. (1997), 
which are based on the two limiting cases of radiative and adi- 
abatic cooling layers behind shock fronts (see above). 

3.6.1. Test calculations 

In the following we present results of test calculations show- 
ing the influence of our modified treatment of shock emission. 
For this purpose we selected the 04f-star C, Puppis as a test ob- 
ject and ignored for the corresponding model calculations the 
improved blocking and blanketing treatment discussed above. 
This restriction makes our results directly comparable to those 
of Pauldrach et al. (1994a), who used the old, simplified treat- 
ment for the shock emission. The stellar parameters of C, Pup- 
pis, used as basic input for our models, have been adopted from 
Pauldrach et al. (1994) (see Table ||) together with the abun- 
dances 

rHc = 1.2OrHc,0 Fc- 0.35 Fee 
Fn^S-OOFn,© Fo- 0.75^0,0 

(Ix '■— nx/riH, where lx,o denotes the solar abundance.) For 
all other abundances solar values were used. 

Although the final objective of our treatment is the deter- 
mination of the maximum post-shock temperature (Tg) and 
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Table 2. The stellar parameters of the 04f-star (" Puppis. 
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Table 3. Parameters required for describing the cooling zones of the 
shocked gas. 
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Fig. 15. Calculated and observed UV spectrum for the 04f-star Puppis. The calculated spectrum belongs to a model where the influence of 
shock emission has been neglected. The high resolution observations have been obtained with the lUE and Copernicus satellites. (Note that the 
improved blocking and blanketing treatment has not been considered for the model calculations of this object — see text.) 
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Fig. 16. Calculated and observed UV spectrum for the 04f-star f Puppis. The calculated spectrum belongs to a model where the influence of 
shock emission has been included — former method (model 1). 
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the filling factor (/) from a comparison of the calculated and 
observed ROSAT spectrum, we have also adopted these val- 
ues for the present test calculations from the similar fits per- 
formed by Feldmeier et al. (1997). The values are given in Ta- 
ble H together with the interstellar column density of hydrogen 
(log(7VH), cf. Shull and van Steenberg, 1985). 

We start with a spectrum synthesis calculation where EUV 
and X-ray radiation by shock heated matter is neglected. The 
comparison between the observed and the synthetic spectrum 
(Fig. |lj) shows clearly that the strong observed resonance lines 
of N V and O VI are not reproduced by the model. This striking 
discrepancy illustrates what is meant by the problem of "su- 
perionization" . In Fig. |l^ we demonstrate, however, that this 
problem has already been solved by making use of the EUV 
and X-ray radiation resulting from the treatment of isothermal 
shocks (previous method — model 1). The observed resonance 
lines of N V and O VI are reproduced quite well, apart from 
minor differences. Thus it seems that the wind physics are cor- 
rectly described. That this is actually not the case can be in- 
ferred from Fig. |l^ where the ROSAT PSPC spectrum (error 
bars) is shown together with the result of model 1 (thin line). 
The deficiency of the non-stratified isothermal shocks is obvi- 
ous — the model yields much too little radiation in the soft 
X-ray part (shortward of 0.7 keV the spectrum is more likely 
characterized by a cooler shock component of logT's ~ 6.30) 
and considerably too much in the harder energy band. 

Following the strategy outlined above we now investigate 
how far the structured cooling zones behind the shocks can 
influence this negative result. Fig. |T7|, which shows in addi- 
tion the calculated X-ray spectrum of our new model (model 2, 
thick line), illustrates the improvement. Strikingly, the new cal- 
culations can quite well reproduce the ROSAT PSPC spectrum 
and the comparison shown is at least of the same quality as 
that obtained by Feldmeier et al. (1997) with their best fit (see 
also Stock, 1998) — note that the total X-ray luminosity of this 
model is given by 



2.5 kcV 



bol- 



(32) 



0.1 kcV 



Actually, it is the fact that, compared to the non-stratified 
isothermal shocks, the post-shock cooling zones with their tem- 
perature stratifications radiate much more efficiently in the soft 
spectral band which leads to the improved fit. This is portrayed 
in Fig. 18 Fig. 19 shows the location of the optical depth unity 
in the relevant energy band of ROSAT. Apart from displaying 
the influence of the K-shell opacities, it becomes evident from 
this figure that the wind is optically thick up to large radii, espe- 
cially in the soft X-ray band. This fact reduces the significance 
of the fit of the ROSAT spectrum, because most of the observed 
X-ray radiation is obviously emitted in the outermost part of 
the wind and thus only the properties of the radiation produced 
in this region can be analyzed from the observed spectrum. 
This, however, is not the case for the EUV and X-ray radia- 
tion which populates the occupation numbers connected with 
the resonance lines of N V and O VI, since due to their P-Cygni 
structure these lines provide information about the complete 
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Fig. 17. Comparison of the ROSAT-observations (error bars) with the 
results of model 2 (thick line) and model 1 (thin line) for Puppis. 




Fig. 18. The Eddington flux in the EUV and X-ray band resulting 
from model 2 (thick line) and model 1 (thin line) for C, Puppis. Note 
that the maximum shock temperatures are identical for both models. 



1000 




E/keV 



Fig. 19. Spatial location of the optical depth unity in the relevant en- 
ergy band of ROSAT, displaying the origin of the observed flux. Also 
shown is the influence of the K-shell opacities for model 2 of ^ Puppis. 
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Fig. 20. Calculated and observed UV spectrum for the 04f-star ( Puppis. The calculated spectrum belongs to a model where the influence of 
shock emission has been included in accordance with our improved method (model 2). 



wind region and the properties of the influencing radiation pro- 
duced in the the whole wind region can therefore be analyzed 
by means of the spectral line diagnostics. Hence, for the signifi- 
cance of our modified method it is therefore extremely convinc- 
ing that the synthetic UV-spectrum resulting from model 2 re- 
produces the observed resonance lines of N V and O VI as well, 
as is shown in Fig. ^ That both model 1 and model 2 yield 
a good fit of the P-Cygni lines shows, on the other hand, that 
distinguishing between two different models from the profiles 
alone is not always possible. The fact that our new treatment 
accounting for the structured cooling zones behind the shocks 
solves not only the problem of " superionization ", but repro- 
duces for the first time consistently the ROS AT PSPC spectrum 
as well as the resonance lines of N V and O VI gives us confi- 
dence in our present approach. (Note that the new treatment of 
the X-ray radiation is not yet available in the download version 
of the code; it will be implemented in an upcoming version 
(2.x).) 

4. Results 

In the following we apply our improved code for expanding 
atmospheres to a basic model grid of 0-stars. The objectives 
of these calculations are to present ionizing fluxes which can 
be used for the quantitative analysis of emission line spec- 
tra of H ll-regions and Planetary Nebulae, and to prove our 
method and demonstrate its reliability by means of synthetic 
UV spectra which are qualitatively compared to correspond- 
ing observations. (Note that for the standard model calcula- 



tions the EUV and X-ray shock radiation is not included — us- 
ing our WM-basic program package this should always be the 
first step. For succeeding models in an advanced stage we have 
used solely our previous method based on isothermal shocks 
(cf. Section 3.6), since this is the method which is presently 
available for WM-basic and thus the models presented in the 
following can be reproduced by this offered tool.) 

Finally, one of the grid models is chosen for a detailed 
comparison between observed and calculated synthetic spec- 
tra, where the primary objective has been to develop diagnostic 
tools for the verification of stellar parameters, and the determi- 
nation of abundances and stellar wind properties entirely from 
the UV spectra. This has been carried out for a cooler 09.5 la 
supergiant, a Cam — a cooler object has been chosen since 
several aspects tend to make these generally more problematic, 
such as the ionization balance (more stages are affected) and 
the optical thickness of the continuum in the wind part. 

4.1. The basic model grid 

In this section we present the ionizing fluxes and synthetic 
spectra of a basic model grid of O-stars of solar metallicity, 
comprising dwarfs and supergiants with effective temperatures 
ranging from 30,000 to 50,000 K. The model parameters, sum- 
marized in Table were chosen in accordance with the range 
of values deduced from observations as tabulated by Puis et 
al. (1996). 

In Fig. |T] we show for each model the primary result, the 
ionizing emergent flux together with the corresponding con- 
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Table 4. The parameters of our basic model grid stars. The Zanstra integrals given here are defined as Qx = / -r— dj^, where hvx is the 

/ hv 

ionization energy of ion X. 



Model 




logs 


R 




M 


logQn 


log (3hc+ 


i^,.(5480A) 
(10^^ erg/s/cm^/Hz) 




(K) 


(cgs) 


(Re) 


(km/s) 


(IQ-^Me/yr) 






Dwarfs 


D-30 


30000 


3.85 


12 


1800 


0.008 


21.42 


8.42 


0.3702 


D-35 


35000 


3.80 


11 


2100 


0.05 


22.65 


11.41 


0.4771 


D-40 


40000 


3.75 


10 


2400 


0.24 


23.15 


17.63 


0.5859 


D-45 


45000 


3.90 


12 


3000 


1.3 


23.45 


18.99 


0.6817 


D-50 


50000 


4.00 


12 


3200 


5.6 


23.69 


20.28 


0.7743 


D-55 


55000 


4.10 


15 


3300 


20 


23.89 


20.17 


0.8881 


Supergiants 


S-30 


30000 


3.00 


27 


1500 


5.0 


22.32 


6.39 


0.4229 


S-35 


35000 


3.30 


21 


1900 


8.0 


22.88 


9.70 


0.4935 


S-40 


40000 


3.60 


19 


2200 


10 


23.19 


11.24 


0.5998 


S-45 


45000 


3.80 


20 


2500 


15 


23.48 


11.84 


0.7160 


S-50 


50000 


3.90 


20 


3200 


24 


23.71 


18.34 


0.8204 



tinuum. (We note that this is the continuum obtained from 
opacities and emissivities resulting from the full line blocking 
and blanketing calculation, which is quite different from the 
one that would be obtained if only continuum opacities were 
used in the iteration cycle.) It can be verified from the figure 
that the influence of the line opacities, i. e., the difference be- 
tween the continuum and the total flux, increases from dwarfs 
to supergiants and from cooler to hotter effective temperatures. 
Both points are not surprising, because they are directly cou- 
pled to the mass loss rate (Af) which increases exactly in the 
same manner (cf. Table Due to the increasing M the op- 
tical depth of the lines also increases in the wind part and 
in consequence the line blocking effect is more pronounced. 
This behaviour, however, saturates for objects with effective 
temperatures larger than Tctf — 45 000 K, since in this case 
higher main ionization stages are encountered (e. g., Fe V and 
Fe VI) which are shown to have less bound-bound transitions 
(cf. Pauldrach, 1987). Thus, as can be verified from Fig. |T] 
the effect of line blocking is strongest for supergiants of inter- 
mediate Tetf . In Table ^ we present the numerical values of the 
integrals of ionizing photons emitted per second for H (log Qh) 
and He II (log (3hc+)> as well as the flux at the reference wave- 
length A = 5480 A, which can be used directly to calculate 
Zanstra ratios and Stromgren radii. 

The next step is to demonstrate the reliability of the cal- 
culated emergent fluxes. As the wavelength region shortwards 
of the Lyman edge usually cannot be observed and thus a di- 
rect comparison of the fluxes with observations is not possible, 
an indirect method to test their accuracy is needed. In princi- 
ple, two such methods exist. The first one is to test the ion- 
izing fluxes by means of their influence on the emission lines 
of gaseous nebulae, i. e., using the ionizing fluxes as input for 
nebular models and comparing the calculated emission line 
strengths to observed ones. However, as a first step this pro- 
cedure is questionable, since the diagnostics of gaseous nebu- 



lae is still not free from uncertainties — dust clumps, complex 
geometric structure, etc. — and therefore, if discrepancies are 
encountered, it is difficult to decide which of the assumptions 
is responsible for the disagreement. (As an example we men- 
tion the Ne Ill-problem discussed comprehensively by Rubin et 
al. (1991) and SeUmaier et al. (1996).) Rather, nebular model- 
ing and diagnostics should be able to build upon the reliability 
of the ionizing fluxes, and thus the quantitative accuracy of the 
fluxes needs to be tested independently of their use in nebular 
emission line analysis. 

The second — and in the light of the difficulties discussed 
above, the only trustworthy — method is quite analogous, but 
instead of an external nebula involves the atmosphere of the 
star itself. The rationale is that the emergent flux is but the outer 
value of a radiation field calculated selfconsistently throughout 
the entire wind, which influences the ionization balance at all 
depths. This ionization balance can be traced reliably through 
the strength and structure of the wind lines formed everywhere 
in the atmosphere. Hence it is a natural and important step to 
test the quality of the ionizing fluxes by virtue of their direct 
product: the UV spectra ofO stars. 

4.2. Qualitative comparison with observations 

The test is performed by means of synthetic UV spectra which 
are qualitatively compared to observed lUE spectra. For this 
comparison we have chosen, for each model of a subset of our 
grid, a real object from the list of Puis et al. (1996) whose sup- 
posed stellar and wind parameters come very close to those of 
the model. The parameters of the model stars and the selected 
real objects are summarized in Table (The influence of shock 
radiation on the models has been neglected at this qualitative 
step.) 

First we investigate the spectra of the dwarf models. As 
can be inspected from Fig. |2| the comparison of the mod- 
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Fig. 21. Calculated ionizing fluxes (Eddington flux in cgs units) versus wavelength of the model grid stars; dwarfs on the left, supergiants on 
the right. 
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Fig. 21. (continued) Calculated ionizing fluxes (Eddington flux in cgs units) versus wavelength of the model grid stars; dwarfs on the left, 
supergiants on the right. 



26 



A.W. A. Pauldrach, T. L. Hoffmann, M. Lennon: NLTE line blocking and blanketing 



D-30 


N IV S VI P IV N IV 

1 


IK Nlll 




Siv Pv Cm Nlll Sim 






^ 








pji if 



Sim Nv Cm 


Oiv Ov Sii\ 









HD 149757 (^Oph) 



Niv SviPivNiv Cm Nlll 



Ovi Siv 



950 1000 1050 1100 1150 1200 1250 1300 900 



950 1000 



1200 1250 1300 1350 1400 1450 1500 1550 



C IV He II N IV 




\imw\i 







Wavelength (A) 





1200 



1050 1100 

Oiv Ov Siw 



1150 1200 1250 1300 




1250 1300 1350 1400 1450 1500 1550 



I 

1500 



1600 1650 1700 1750 

Wavelength (A) 



1800 1850 



' — 
900 



D-40 



' S VI P IV Niv Chi Nim Ovi Siv 



Cm Nlll Sim Nv Cm 




950 1000 



1050 1100 

O IV O V Si w 



1150 1200 1250 1300 




HD 217068 



Niv SviPivNiv Cm Nm 



Ovi Siv 





900 950 




1050 1100 

O IV O V Si w 



1150 1200 1250 1300 



1200 1250 1300 1350 1400 1450 1500 1550 1600 1200 1250 1300 1350 1400 1450 1500 1550 1600 



1500 1550 



1650 1700 1750 

Wavelength (A) 



1800 1850 1900 1500 1550 1600 



1650 1700 1750 

Wavelength (A) 



1800 1850 1900 



D-50 





vj It IV Niv Cm Niii Ovi Siv Pv Cm Niii Sim N 

1 \ h ^ . . . . 


K „■ 











900 
2 



950 1000 1050 1100 1150 1200 



1250 1300 



Si III N 


V C III IV C 


V Si IV Sv C 


K - 


\ J 


l^^*^rt¥^l^t,"y|,^l\l\| 


, i/ 





1200 1250 



1300 1350 1400 1450 1500 1550 



HD 93250 



Niv SviPivNiv Cm Niii Ovi Siv Pv 













1000 1050 1100 1150 1200 1250 



Si Nv C 


III Oiv Ov Siiv Sv Civ 

, . 1 








u 


mm 



1600 1200 1250 1300 



1350 1400 



1450 1500 1550 1600 



\\.ffl'. ' . ! J I 



1500 1550 



1650 1700 1750 

Wavelength (A) 



1800 1850 1900 1500 1550 1600 



1650 1700 1750 

Wavelength (A) 



1800 1850 



Fig. 22. Calculated UV spectra of the model grid dwarfs (left) compared to observed lUE spectra (right) of stars of similar spectral type. 
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Table 5. Model grid stars and real-world examples with similar spectral types. The parameters of the observed example stars are from Puis et 
al. (1996). 



Model / 




logs 


R 




M 


Example 


(K) 


(cgs) 


(Rq) 


(km/s) 


(W-H4Q/yr) 


Dwarfs 


D-30 


30000 


3.85 


12 


1800 


0.008 


HD 149757 (( Oph) 


32500 


3.85 


12.9 


1550 


<0.03 


D-40 


40000 


3.75 


10 


2400 


0.24 


HD 217068 


40000 


3.75 


10.3 


2550 


<0.2 


D-50 


50000 


4.00 


12 


3200 


5.6 


HD 93250 


50500 


4.00 


18 


3250 


4.9 


Supergiants 


S-30 


30000 


3.00 


27 


1500 


5.0 


HD 30614 (a Cam) 


30000 


3.00 


29 


1550 


5.2 


S-40 


40000 


3.60 


19 


2200 


10 


HD 66811 (C Pup) 


42000 


3.60 


19 


2250 


5.9 


S-50 


50000 


3.90 


20 


3200 


24 


HD93129A 


50500 


3.95 


20 


3200 


22 



els D-30 and D-40 with their counterparts HD 149757 and 
HD 217068 show in principle an overall agreement, whereas 
the D-50 model, compared with its counterpart HD 93250, 
shows a severe discrepancy concerning the Ov subordinate 
line at 1371 A (the calculated line is much too strong) and a less 
pronounced discrepancy of the N IV subordinate line at 1718 A 
(the calculated line is somewhat too weak). Hence we have to 
realize that either the wind physics is not completely described, 
or the stellar or the wind parameters of this model are too dif- 
ferent from those of HD 93250. 

Regarding the first point one might speculate that the in- 
clusion of shock radiation leads to an improvement for the O V 
line, although this effect would weaken the N IV line further 
As is shown below, shock radiation cannot solve the problem, 
as it does not affect the strength of the O V line at all (cf. the 
discussion of the S-50 model below). Regarding the second 
point there are three parameters which could lead to an im- 
provement for both lines. The first one is the effective temper- 
ature which, however, would have to be decreased by at least 
5000 K. This is on the one hand extremely unrealistic, since 03 
and 04 stars would have almost the same Tcff, and it would on 
the other hand produce another discrepancy due to an increase 
of the strength of the Oiv fine at 1338 A (cf. the S-40 model 
in Fig. |4|). The second parameter is the mass loss rate and 
the third one is the abundance. In order to investigate whether a 
systematic variation in the mass loss rate can solve the problem 
we computed a small model subgrid for this object by changing 
the mass loss rate for model D-50, keeping all other parameters 
the same (cf. Table 

The synthetic spectra obtained for the mass loss rate se- 
quence are shown in Figure As can be seen, lowering the 
mass loss rate does not solve the problem of the too-strong 
O V line, as this is one of the last lines to disappear with di- 
minishing mass loss rate, whereas the N IV line which was al- 



ready too weak disappears immediately. A mass loss rate as 
low as 10~* Mq/yi would be required to reduce the strength 
of the O V line to the observed case. The only model where 
the situation has clearly improved regarding both the N IV and 
the O V line is D-50-a, the one with a higher mass loss rate. 
Hence we conclude that the mass loss rate has to be larger by 
a factor of « 2 — note that this is also indicated by the stellar 
parameters which are almost identical to those of HD 93 129 A, 
which means that the wind parameters have to be almost identi- 
cal too; note further the strong similarity of the observed spec- 
tra of HD 93 129 A and HD 93250, which points in the same 
direction (cf. Fig. ^s], upper panel); in addition, the abundances 
of the CNO elements have to be reduced strongly compared to 
the solar values (for oxygen a reduced value of approximately 
a factor of 50 is required to weaken the strength of both the 
O IV and the O V line, as has been inferred from additional 
test calculations). This, however, should not be the case for 
the heavier elements (Fe, Ni) since a solar-like abundance of 
these elements is required to account for the radiative acceler- 
ation necessary to produce the higher mass loss rate (cf. Paul- 
drach, 1987). 

Now we inspect the comparison of the supergiant models 
(Fig. |4|). As is shown, the observed spectra are reproduced in 
principle quite well apart from minor differences which can be 
attributed to a change of abundances (note that the discrepancy 
of the N V resonance line is due to the omitted shock radiation 



Model 


D-50-a 


D-50 


D-50-b 


D-50-C 


D-50-d 


M(lO-«M0/yr) 


11.0 


5.6 


0.56 


0.12 


0.01 



Table 6. Mass loss rates for the models shown in Figure g^. All other 
parameters are identical to those of the D-50 model shown in Table H. 
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Fig. 23. Spectra obtained for a mass loss rate sequence for the D-50 dwarf model. For comparison, the lUE spectrum of HD 93250 is also 
shown. The mass loss rates of the models are given in Table ^. Lowering the mass loss rate cannot solve the problem of the too-strong O V line, 
as this is one of the last lines to disappear with diminishing mass loss rate. 
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Fig. 24. Calculated UV spectra of the model grid supergiants (left) compared to observed lUE spectra (right) of stars of similar spectral type. 



30 



A.W. A. Pauldrach, T. L. Hoffmann, M. Lennon: NLTE line blocking and blanketing 



HD 93129A&HD 93250 



HD93129A 
HD 93250 



N IV S VI P IV N IV Cm N III 



1.5 

1 

0.5 




900 950 1000 1050 1100 1150 1200 1250 1300 



2 

1.5 

1 1 
ai 

0.5 



Nv.p 



IV V Si IV 



^ i|,i, '1 ■ ' - ■ 




1900 



900 950 1000 1050 1100 1150 1200 1250 1300 



2 

1.5 - 
0.5 - 



Sim Nv Cm 


IV V Si IV 


Sv C 















1200 1250 1300 1350 1400 1450 

2 I n 1 , 

1.5 - 



1500 1550 1600 



1 p'VT ^ pr""-' YvYT" 



I — 
1500 



1550 1600 



1650 1700 
S-50-a 



1750 1800 1850 1900 



1.5 - 



Niv SV|R 



C III N III SUM N V C III 



^ If 

0.5 - 




900 950 1000 1050 1100 1150 1200 1250 1300 



2 
1.5 

1 

0.5 




Si HI 


N 


V cm Oiv 


Ov Siiv Sv C 




- 


/ 











1200 1250 1300 1350 1400 1450 

2 I n 1 1 



1500 1550 1600 



1.5 - 
1 



0.5 - 



.-.A— V — w /N. 



1500 1550 1600 1650 1700 1750 1800 1850 1900 

Wavelength (A) 

Fig. 25. Upper panel: Comparison of the lUE spectra of HD 93129A (03 If) and HD 93250 (03 V ((f))). Apart from the strength of a few 
lines of light elements (e. g., the N IV line at 1718 A) the spectra are almost identical. Other panels: Influence of shocks on the spectrum of the 
S-50 supergiant model which is compared to the lUE spectrum of HD 93129A. Middle panel: without shocks, lower panel: with shocks. The 
influence of the shocks on the strength of the O IV and O V lines is negligible. 




Table 7. Shock parameters for the S-50-a model — for an explanation 
of the parameters see Section 4.3. 

log(Lx/I/boi) Vt/voc, 7 m 

^7i) oi i i~ 



— see below). Again, the most conspicuous difference regards 
the O IV and O V subordinate lines which are both too strong, 
especially for the S-40 and the S-50 models. From the inves- 
tigation above it is akeady quite clear that the abundance(s) 
of the (CN)0 element(s) has (have) to be reduced in order to 
overcome this discrepancy. Nevertheless, we investigate now 
whether the inclusion of shock radiation leads to an improve- 
ment. We have therefore computed an additional model for the 
S-50 supergiant, model S-50-a, where the influence of shocks 
on the spectrum has been accounted for; the shock parameters 
are given in Table ^ 

Fig. ^ shows that the influence of the shocks on the 
strength of the O IV and O V lines, and hence the ionization bal- 
ance, is negligible — just O VI is enhanced selectively, which 
can be seen by the strength of the Ovi resonance line. An- 
other line which is considerably affected by shock emission is 
the Si IV resonance line. This is because the soft X-ray radi- 
ation field of the shocks enhances the ionization of Si V, and 
thus the recombination to Si IV is decreased (cf. Pauldrach et 
al., 1994). As can be inferred from Fig. ^ this improves the fit 
of the Si IV line significantly. Concerning oxygen we come to 
the same conclusion as for the dwarf models, namely that the 



abundance has to be reduced strongly compared to the solar 
value. 

4.3. Detailed analysis of a Cam 

In this section we provide, using our S-30 grid model as a start- 
ing point, a detailed determination of the abundances and stel- 
lar wind properties, and verification of the stellar parameters 
of a Cam. Special emphasis is given to the shock radiation 
needed for a fit of the N V and O VI resonance lines. Although 
preliminary results fro m ou r new shock description look very 
promising (see section |3.6| ), this new method is not yet fully 
implemented in WM-basic. For this investigation we therefore 
use the method based on isothermal shocks. See Pauldrach et 
al. (1994) for a detailed explanation of the shock parameters 
and the rationale behind the parametrization. 

We wish to point out here that we have not attempted an 
exact determination of log g and stellar radius, but have rather 
kept the values of our S-30 grid model. The reason is that in 
contrast to the hydrodynamics the UV spectrum depends only 
marginally on these parameters, the main influence being due 
to Tcff and M (i. e., density). The radius and the surface gravity 
(log g) can in principle be determined from a selfconsistent cal- 
culation of the hydrodynamics and the NLTE model, in which 
both values would be adapted in such a way that the hydro- 
dynamics, with consistent force multiplier parameters from the 
NLTE occupation numbers (in turn again dependent on the hy- 
drodynamic solution), would yield the mass loss rate and the 
terminal velocity deduced from Ha and the observed UV spec- 
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Table 8. Model parameters for a Cam. Abundances not explicitly mentioned are solar; if given, the number is the factor relative to the solar 
value. Toff is in K, M in 10"*^ Mg/yr. 



Model 


Tcff 


M 




log(I/x/ibol) 


7 


m 




abundances 




a 


30000 


5 


0.125 


-7.5 


1 


1 




solar 




bl 


30000 


5 


0.14 


-7.5 


1 


1 




solar 




b2 


30000 


5 


0.14 


-6.0 


1 


1 




solar 




b3 


30000 


5 


0.14 


-6.5 


1 


1 




solar 




c 


30000 


5 


0.20 
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0.5 




C = 0.1,N 


= 10., 0= 1.0, P 


= 0.1, S = 2.0 


dl 
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' 11 
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h 


28000 


5 


0.25 


-6.5 


0.5 


30 
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trum (cf. Pauldrach et al., 1994, 1994a). This latter procedure 
has however not yet been implemented in WM-basic. 

In Fig. ^ we compare the spectrum of our basic S-30 grid 
model with the spectra of a Cam observed with lUE and Coper- 
nicus. With the exception of a few strong lines, notably the N V 
resonance doublet at 1238, 1242 A and the subordinate Cm 
line at 1247 A, the agreement is very good for a first step, espe- 
cially with respect to the iron and nickel "forest" between 1400 
and 1600 A. As the two lines mentioned above are mainly af- 
fected by shock radiation, we will first attempt to fit the shock 
parameters. We note in passing that the strong lines from 1000 
to 1100 A (with the recurring pattern), that complicate an ex- 
act fit of the O VI resonance line, are absorption by interstellar 
molecular hydrogen. 

A model with reasonable "first guess" values for the shock 
parameters (model a. Table ^ already gives a very good fit to 
the N V resonance line, as shown in Fig. the emission, how- 
ever, completely decimated by the much too strong C III A 1247 
line. Increasing the X-ray luminosity {log{Lx/ L]^o\)) to reduce 
the C III occupation unfortunately also tends to ionize Si III 
to Si IV, thus reducing the strength of the Si III line. Modi- 
fying vt /voo (the ratio of the maximum jump velocity to the 
terminal velocity, where t^t characterizes the immediate post- 
shock temperature; see Pauldrach et al., 1994) does not change 
this fact, as we have confirmed by calculating a grid of mod- 
els with values of log(Lx/iboi) ranging from —8 to —6.5 
and Vt/voo ranging from 0.1 to 0.2. For example, a model se- 
quence (models bl-b3) is shown in Fig. ^ in which we vary 
log(Xx/iboi) from —7.5 to —6.5, keeping all other parameters 
constant (cf. Table |^). As can bee seen, the Si III line already 
begins to weaken, while C III still remains too strong. Adjust- 
ing the 7-parameter (which controls the strength of the shocks 
relative to the local velocity — cf. Pauldrach et al., 1994) does 
not help either, as both the C III and the Si III line are formed in 



the inner part of the wind, thus both being subject to the same 
radiation field. 

From this we conclude that solar abundances cannot repro- 
duce the observed spectrum. A model, however, with a reduced 
carbon abundance of one tenth solar (model c) can indeed give 
a good fit of the line, as shown in Fig. ^ We will use this 
abundance for the following calculations, unless stated other- 
wise. As carbon thus shows indications of the CNO-process, 
we have for reasons of consistency also increased the nitrogen 
abundance by a factor of 10, but such a large factor is not com- 
patible with the fit of the nitrogen lines; from the final models 
we determined the nitrogen abundance to be approximately so- 
lar. 

Despite the good fit, this model is still not satisfactory, as it 
shows no signs of O VI. Since the C Ill/Si Ill-balance strongly 
constrains the shock strength in the inner regions, this cannot 
be achieved simply by increasing the X-ray luminosity. How- 
ever, the onset of shocks can be adapted with the m-parameter, 
which gives the ratio of outflow to sound velocity where shocks 
start to form, and as the jump velocity is correlated to the 
outflow velocity, the corresponding radius (cf. Pauldrach et 
al., 1994). Increasing this parameter allows shocks in the outer 
regions (Ovi appears close to Voo), while leaving the inner 
regions (where Si III is present) largely undisturbed. This be- 
haviour of the shocks is already an important result from our 
analysis. 

Fig. ^ shows a model sequence in which m is increased 
from 1 to 60; \og{Lx/ L}^o\) has been increased compared to 
model c to provide for sufficient O VI. It can bee seen that in 
model d4 (m = 60), where Si III retains its correct strength and 
the absorption of O VI and N V is of the correct magnitude, the 
emission of both latter lines has disappeared completely. This 
is because the emission of P Cygni lines due to resonance scat- 
tering arises to a large part from the inner regions of the wind, 
where the stellar radiation field is still strong; in this model 
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Fig. 27. Model a, using our "first guess" values for the shock parameters. C III A 1247 is much too strong. 
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Fig. 28. Models bl-b3, a sequence showing the influence of increasing log(Lx/I'boi). Note that Sim disappears before C III has decreased 
to its observed strength. 
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Fig. 29. Model c. Here we have adapted the shock parameters and abundances (see Table H) to provide for a reasonable fit of Si III and C III. 
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Fig. 30. Models dl-d4, a sequence in which the radius at which shocks start to form has been successively moved outwards, keeping all 
other parameters equal. In model d4, where Si ill remains fairly unaffected by shocks, the P Cygni emission of O VI and N V has completely 
disappeared. 
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Fig. 31. Model e. Increasing the shock jump velocity provides for harder shock radiation that has a stronger influence on O VI and N V than on 
Si III. The abundances of C, N, O, and P have also been adapted. 
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Fig. 32. The complete EUV spectrum of model e (Teff = 30000 K). 
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Fig. 33. Model f. Like model e, but with a Teff of 29000 K. Note the better fit of Si ill and O IV. 
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Fig. 34. Model g. Like model f, but with adapted abundances (see Table B). This is our best fit. 
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Fig. 35. Model h. Same as model g, but with T^fc = 28000 K 
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Fig. 36. Model i. As model h, but with a lower mass loss rate (M = 2.5 x 10 ® M© /yr compared to M = 5 x 10 Mq /yr for model h). 
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Fig. 37. Model j. As before, but with a higher mass loss rate (M = 10 x 10 ® Mq /yr). Note the change in the strength of the He ii line. 



38 



A.W. A. Pauldrach, T. L. Hoffmann, M. Lennon: NLTE line blocking and blanketing 



S-30 + model-e 



S-30 
model e 



E 



S-10 

LL 

o 

■5)-12 
g 

T3 
LU 



-14 
-16 
-18 

-20 



100 200 300 400 500 600 

Wavelength (A) 









































r^ir^nf 








.vipilf 


ilf* 


jRfli iff f T 

ll ' ^l' 1 


1 


11 
ll 


1 1 f 
1 1 " 








i||f1|t, 

l(| 

1'' 


10: 




1 ' i ' 




1 


























III ll 








■if 


f 


















I 




















1 







































































700 800 900 1000 



Fig. 38. Comparison of the ionizing flux of our S-30 grid star, wliere shocks have been neglected, to that of model e, showing the non-negligible 
influence of shocks on the ionizing tiuxes. 



however, N V and O VI hardly exist in these regions, since the 
shock radiation responsible for the presence of these ions is 
only produced far out in the wind and the ionization contin- 
uum is still optically thick at the corresponding photoioniza- 
tion thresholds. We conclude that the onset of shocks must lie 
further in, and, correspondingly, the shock radiation must be 
harder (maximum at smaller wavelengths) than attainable with 
a Vt/voo value of 0.14, so as to have a comparatively larger in- 
fluence on N V and O VI via Auger-ionization than on Si III via 
photoionization. 

This reasoning is confirmed by model e (Fig. in which 
we have reduced m to 30 and increased Ut/^^oo to 0.25 (cor- 
responding to a maximum shock temperature of 2.0 x 10^ K). 
The shock parameters of this model yield the best overall agree- 
ment of the most important UV spectral lines. Note that for this 
model we have in addition reduced the phosphorus abundance 
to 0.05 times solar, to obtain a fit of the Pv line. As it con- 
cerns only a single element and since evidence for an under- 
abundance of phosphorus has already been encountered from 
an analysis of another hot star, namely the 04f-star ( Puppis 
(cf. Pauldrach et al., 1994; and Section 3.6, this paper), it is 



most likely that the discrepancy between the observed and cal- 
culated P V resonance line is caused by the proposed abundance 
effect. We further note that the imbalance of S IV and S V is 
most likely due to a combination of abundances and imper- 
fect atomic data, as our S IV atomic model has not yet reached 
the quality of those of other ionization stages and is still in a 



stage of rather incomplete description (cf. Table Note that 
we have also reduced the oxygen abundance, thereby improv- 
ing the fit of both O IV and O VI. 

Having thus constrained the strength and distribution of 
the shocks, we still need to check to what extent the effective 
temperature and the mass loss rate can be constrained further 
through our analysis of the UV spectrum. For this purpose, we 
have computed a model with Tcf f of 29 000 K but otherwise 
same parameters (model f), whose spectrum is shown in Fig- 
ure [33I C omparing this to the spectrum of model e (30 000 K) 
(Fig. [32[) we note a marginally better fit of the spectral region 
from 1450 to 1650 A, and a slight improvement in Si IV and 
Oiv. More significant is the extreme sensitivity of the Sim 
and the Cm line to this small change in temperature, which 
reveals that these lines can be utilized as temperature indica- 
tors in this spectral range. We concede that the fit of N V is still 
somewhat imperfect, but point out that the shock model used is 
not in a final stage (see above); we expect this to improve with 
our new method. In a next step we have reduced the carbon 
abundance further (to 0.05 solar) and nitrogen to solar abun- 
dance, as this improved not only the N v/C III fit, but also the 
N III fine at 1 183,1 185 A (Fig. |^). Note also that the saturated 
C IV resonance line is not affected by this change of the carbon 
abundance. 

Lowering the effective temperature further to 28 000 K 
worsens the fit, as the iron and nickel lines around 1500 A be- 
come too strong (Fig. ^s]). As the strength of these lines de- 
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pends also on the mass loss rate, it is conceivable that an adap- 
tation of this parameter can again improve the fit. However, 
M is strongly constrained by the strength of the Hell line, 
as demonstrated in figures ^ {M — 2.5) and|3^ {M ~ 10). 
We have confirmed in test calculations that the iron abundance 
must be solar {Z = 1.0) to reproduce the observed spectrum 
and to account for the radiative acceleration needed to produce 
the observed mass loss rate (cf. Pauldrach, 1987), a finding 
compatible with our earlier statement concerning the relative 
abundances of iron to the light CNO elements. 

Finally, we wish to demonstrate the influence of the shocks 
on the ionizing fluxes. In Figure 38 we compare the ionizing 
flux of our S-30 supergiant grid model, where shocks have been 
neglected, to that of our model e for a Cam. (We have used 
model e for this comparison, and not our final model g, since 
model e and the grid star S-30 both have the same effective 
temperature of 30 000 K.) Due to the lower optical depth red- 
wards of the He II edge, the shock radiation enhances the flux as 
far as redwards as 400 A; this will have important implications 
concerning the solution of the Ne Ill-problem (cf. Sellmaier et 
al. 1996), pending further investigation. Note: the Zanstra inte- 
gral for He II is increased from log (5ho+ = 6.39 for the S-30 
model to 15.07 for model e! 



5. Conclusions 

After a long period of work in the areas of non-LTE radia- 
tive transfer, hydrodynamics, and atomic physics we have now 
developed a fast numerical model code for expanding atmo- 
spheres which incorporates for the first time the physics re- 
quired — rate equations for all ions using detailed atomic 
physics, the equations of radiation hydrodynamics, the energy 
equation, the radiative transfer equation including the effects 
of overlap of numerous spectral lines of different ions, and 
shock emission of instabilities in the stellar wind flow. We have 
shown that modelling hot star atmospheres involves replicating 
the tightly interwoven mesh of these physical processes, the 
most complicating effect in this system being the overlap of 
thousands of spectral lines of different ions. Especially con- 
cerning this latter point we have made significant progress; the 
decisive factor has been to relax some rather severe approx- 
imations concerning the correct treatment of Doppler-shifted 
line radiation transport and the corresponding coupling with 
the radiative rates in the rate equations. We have demonstrated 
that these modifications to the models lead to changes in the 
energy distributions, ionizing continua, and line spectra with 
much better agreement with the observed spectra. This has im- 
portant repercussions for the quantitative analysis of hot star 
spectra. With this new method in hand we have already estab- 
lished a new basic model grid of stars of solar metallicity that 
can be used as input for the analysis of spectra of emission line 
nebulae. The qualitative diagnostic investigation performed on 
basis of the model grid revealed further that the oxygen abun- 
dance is considerably reduced compared to the solar value, at 
least for objects younger than 05. 



From the first detailed analysis of the 09.5 supergiant 
a Cam we conclude that our spectrum synthesis technique 
does, in principle, allow the determination of effective temper- 
ature and abundances — in fact, a determination of the effec- 
tive temperature to within ±1000 K and of the abundances to 
within a factor of 2 seems not unreasonable. Carbon and phos- 
phorus show clear signs of an underabundance on the order of 
one tenth its solar value, as does oxygen with about 0.3 so- 
lar, whereas the abundance of iron must be roughly solar to 
reproduce the spectrum of the numerous Fe IV and Fe V lines. 
To produce the ionization balance observed in the lighter ele- 
ments C, N, O, and Si, the influence of shock radiation must 
start at larger radii where for shorter wavelengths the largest 
shock temperatures dominate. Thus, the way the X-ray spec- 
tral region selectively affects the ionization balance of differ- 
ent elements, observable through the lines in the EUV spec- 
trum, provides constraints on the lower shock temperatures; we 
have determined maximum shock temperatures on the order of 
2.0 X 10^ K. Especially the Si III and C III fines have been found 
to be invaluable diagnostic instruments for this purpose. 

Our detailed analysis of the UV spectrum and the shocks 
needed to reproduce the observed lines has led to a significant 
difference in the ionizing flux of this model. Thus we conclude 
that this type of analysis is indispensable and must be regarded 
as the ultimate test for the accuracy of ionizing fluxes from 
models. 

Our research plan for the future has three major objec- 
tives. First, we will have to implement further improvements to 
the model atmosphere code, especially concerning the planned 
analyses of optical lines. For this purpose Stark-broadening has 
to be included for the affected spectral lines (e. g., H and He I 
lines) and concerning the rate equations instead of using the 
Sobolev-plus-continuum method some of these lines should be 
treated in the comoving frame, if they are used for diagnostic 
purposes (cf. Sellmaier et al., 1993, who found that compared 
to the comoving frame treatment the Sobolev-plus-continuum 
approximation leads to a non-negligible change of the strength 
of some H and He lines). In connection with this some minor 
approximations will also have to be checked in detail. 

Second, we plan to apply the model atmosphere code to 
a comprehensive sample of stars of different metallicities for 
modelling ionizing fluxes and line spectra, the qualitative re- 
view of the model grid presented in this paper being just a first 
step for a detailed quantitative analysis of ionizing fluxes and 
quantitative verification of the accuracy and reliability of the 
models. Furthermore, spectral analysis in the UV and the op- 
tical range of individual stellar objects in Local Group galax- 
ies and in galaxies as distant as the Virgo Cluster will be per- 
formed. 

Third, with the new diagnostic tool of detailed expanding 
model atmospheres presented, and the observations of the HST 
space observatory and the ESO VLT ground-based telescope 
which are already available, the concept of using luminous hot 
stars for quantitative UV spectral analyses for determining the 
properties of young populations in galaxies is not only reason- 
able, but first steps in this direction have already been taken 
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(see Mehlert et al., 2000, who found that the spectra of galax- 
ies they had observed at high redshifts (z ~ 3) displayed the 
typical features usually found in the UV spectra of hot stars; 
first diagnostic investigations have shown that these spectra, 
considering proper reddening, can be well fitted with synthetic 
spectra from our hot star models). Thus, the calculation of inte- 
grated spectra for the determination of stellar abundances and 
the physical properties of the most UV luminous stars in star 
forming galaxies even at high redshifts using the technique of 
population synthesis is feasible. 
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Appendix A: A new concept for a fast solution of the 
Rybicki-metliod 

Here we present for the solution of the Rybicki-scheme a concept 
which is compared to the standard procedure 10 times faster on a 
vector processor and 3 to 5 times faster on a scalar processor (see 
Sec. 3.3). 

We start from the final system for the solution of the mean inten- 
sity — the vector J describes its depth variation (the number of depth 
points is A''): 

J = W'^Q 

with (TV' is the number of p-rays) 

at' 



The usual, but time consuming solution method is to calculate the in- 



W = W 
where 



Q 



n' 



(A. I) 



(A.2) 



rr' -1 



V 



(A.3) 



and Nj is the number of radius points for the jth p-ray (cf. Fig. 
Wj is the diagonal matrix of the integration weights, and Tj is a tri- 
diagonal matrix defined by the coefficients of the difference equation 
of transfer 

UijUi-ij + bijUi^j + CijUi+i^j — PiJi = (1 — Pi)Si (A.4) 

where ai.j, bij, and Cij are the coefficients as in eq. ^ (section 3.3.1) 
for the jth p-ray, 

Cl,.7 \ 

(A.5) 



T' 

^3 



The variables Uij are the symmetric averages of the intensities, the 
coefficients f3i are the ratios of Thomson-opacities to total opacities, 
and the diagonal matrix Uj and the vector Kj are defined as 

/-/9i \ / (l-/9i )Si 



Uj = 



(A.6) 



verse matrices T. 



- by a forward-elimination and back-substitution 



procedure — and from these, obtaining Q and W, and, finally, J. 

Due to the diagonal character of the matrix Uj this is, however, 
not necessary, because from the solution of the first column of W 
obtained by using just the first column of Uj in eq. k.2 the solutions of 



the remaining columns of W can be generated. Hence, the elimination 
procedure only has to be applied to the following two sets of equations 
{N[ is the maximum number of p-rays for the ith radius point): 
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K, -- 
where 
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(A.7) 



(A.8) 



(A.9) 



(A. 10) 



Note that since the solution procedure is not recursive with respect to 
the index j the elimination can be performed simultane ously fo r all 
p-rays. Note further that the structure of the sums in eqs. A.7 and A.J 
are equivalent to the operations x = A ■ y where A is a triangular 
matrix. Thus, BLAS IeveI-2 routines can be applied. 

For the construction of the remaining columns of W we now make 
use of the already calculated matrix Bj^i and two auxiliary matrices 
obtained as a by product during the forward-elimination procedure by 



solving eq. \.8 



The first one is: 
with 

-Pj.O = 1, fj,i = 



i = 1, 



-a 



'j,i + l 



(A. II) 



(A. 1 2) 



Here bjj is the updated valu e of bj,i obtained after the forward- 
elimination step by solving eq. A.8 (&j,i+i = 6j,i+i + fj,iCj,i)- 
The second one is: 



with 



F: 



J,N 



F", 

^ N'.,i 



i = N -I, 



93,: 



'^3,ilbj 



(A. 1 3) 



(A. 14) 



The final step — construction of the columns I = 2, . . . , A'^ of 
W (w here I denotes the position of the unity value for the vector Bj 
in eq. k.li ) — consists of two parts, where the first part replaces the 
forward-elimination and the second part the back-substitution proce- 
dure: 

In the first part the components i = /,..., TV of the l\h column of 
W are determined by 



N'. 



(A. 1 5) 



j=i 
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where 



(A. 16) 



Again, a triangular matrix has simply to be multiplied by a vector 
using a BLAS level-2 routine. 

In the second part the components i — I — 1, ... ,1 of the Zth 
column of W are determined by 



where 

v-:i^vi,b,,/f;:,. 



(A. 17) 



(A. 18) 



Now, a rectangular matrix has simply to be multiplied by a vector 
again using a BLAS level-2 routine. Hence, instead of a number of op- 
erations ~ A'''^ a number of operations only ~ A'^^ must be performed, 
and as extremely fast routines are used to solve the non-recursive sys- 
tem, the solution of the Rybicki-scheme is now almost as fast as the 
solution of the moments equation. 



Appendix B: A simple demonstration of the failure of the 
standard r discretization of the equation of 
transfer 

Let us assume a one-dimensional, sharply bounded cloud of material 
of moderate opacity. 



1, 

10" 



< 5 
> 5 



(B.l) 



with a symmetric boundary condition at z = (the middle of the 
cloud) and no influx from the outside (i. e., = at 2 = 10). 

With a constant source function S = 1 everywhere, both the old 
and the new discretization of the transfer equ ation yield essentially 
the same radiation field (upper panel of Figure |B.1[ ), a solution which 
is immediately obvious: The center of the cloud (z = 0) is optically 
thick, the intensity there thus equal to the source function. Towards the 
edge of the cloud, more and more radiation escapes, and the intensity 
decreases. On the outside of the cloud, essentially no emission is pro- 
duced (nothing is absorbed, either); only the radiation emitted by the 
cloud contributes to the local intensity (since we have no influx from 
the right), therefore u — ^S. 

The results change considerably if we adopt a high source func- 
tion of 5* = 100 for the single grid point at z = 5, the first grid point 



outside of the cloud. As illustrated in the lower panel of Figure B.l 
the old method produces extra emission on a scale of 



AI = -S- At 
2 



ilOO-i(l + 10-"')Az 

ilOOAz. 

4 



For the grid with 100 points, Az = 0.1, giving us an intensity 

11 f» 0.5 + i ■ ilOO-0.1 = 1.75 

2 4 



(B.2) 
(B.3) 

(B.4) 



outside the cloud. (0.5 is from the normal emission of the cloud as 
before, and the factor 1 /2 in the second term is again due to the fact 
that u = ^{I^ + with /~ = on the outside of the cloud, so 
that u — §/"''.) The extra emission is proportionally reduced in the 
finer grids as Az gets smaller. 




Spatial coordinate z 





old method, 100 points 




old method, 200 points 




old method, 500 points 




old method, 1500 points 








new method {all curves) 



Spatial coordinate z 

Fig. B.l. Upper panel: Mean intensity in a one-dimensional nebula; 
constant source function, high opacity on the left (shaded in gray), 
low opacity on the right. Both discretization methods give the same 
results. Lower panel: If the first grid point with low opacity is given a 
high source function, the old method produces artificial emission. The 
new method does not show this behavior. 



The new method, on the other hand, is essentially indifferent to 
the grid spacing in this particular configuration, since the emission is 
only computed on the basis of the local opacity, which is low for the 
point in question. If we were to choose a point inside the cloud to 
have a high source function, then both methods again produce similar 
results, also dependent on grid spacing, since in this case the point in 
question does have a high opacity, and therefore produces substantial 
emission proportional to the interval length. 

Concluding, it can be said that the discretization methods differ in 
the assumptions that are made about the emissivity (in the old method, 
based on average values of opacity and source function assuming lin- 
ear run between the grid points; in the new method, extrapolating the 
local value) — information not present given just the values at the 
grid points, and which must therefore be supplied separately through 
the choice of discretization coefficients. It just so happens that the new 
method produces results which are more compatible with our detailed 
formal integral. 
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